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ABSTRACT 


There  are  numerous  numerical  methods  for  solving  different  types  of  partial  dif¬ 
ferential  equations  (PDEs)  that  describe  the  physical  dynamics  of  the  world.  For 
instance,  PDEs  are  used  to  understand  fluid  flow  for  aerodynamics,  wave  dynam¬ 
ics  for  seismic  exploration,  and  orbital  mechanics.  The  goal  of  these  numerical 
methods  is  to  approximate  the  solution  to  a  continuous  PDE  with  an  accurate  dis¬ 
crete  representation.  The  focus  of  this  thesis  is  to  explore  a  new  Discontinuous 
Galerkin  (DG)  method  for  approximating  the  second  order  wave  equation  in  com¬ 
plex  geometries  with  curved  elements.  We  begin  by  briefly  highlighting  some  of 
the  numerical  methods  used  to  solve  PDEs  and  discuss  the  necessary  concepts  to 
understand  DG  methods.  These  concepts  are  used  to  develop  a  one-  and  two- 
dimensional  DG  method  with  an  upwind  flux,  boundary  conditions,  and  curved 
elements.  We  demonstrate  convergence  numerically  and  prove  discrete  stability  of 
the  method  through  an  energy  analysis. 
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CHAPTER  1 : 

An  Introduction  to  Solving  Partial  Differential 

Equations 


There  are  numerous  numerical  methods  for  solving  different  types  of  partial  dif¬ 
ferential  equations  (PDEs)  that  describe  the  physical  dynamics  of  the  world.  The 
goal  of  these  numerical  methods  is  to  approximate  the  solution  to  the  continuous 
PDE  with  a  discrete  representation  [1],  Three  notable  methods  are  Finite  Difference 
Methods  (FDMs),  Finite  Volume  Methods  (FVMs),  and  Finite  Element  Methods 
(FEMs).  Two  common  FEMs  are  Continuous  Galerkin  (CG)  and  Discontinuous 
Galerkin  (DG),  each  of  which  comprises  an  element-based  approach  to  solving  a 
set  of  equations.  The  main  focus  of  this  thesis  is  to  explore  a  new  DG  method,  in¬ 
troduced  by  Appelo  and  Hagstrom  [2],  for  approximating  the  second-order  wave 
equation.  However,  my  research  differs  from  that  of  Appelo  and  Hagstrom  by 
using  a  nodal  form  of  DG  with  provable  stability  on  curved  elements  with  complex 
geometries. 

1.1  Finite  Difference  Method 

Finite  difference  methods  are  one  of  the  simplest  and  oldest  methods  for  solving 
partial  differential  equations  [1].  Furthermore,  it  is  arguably  one  of  the  most  used 
methods  for  discretizing  partial  differential  equations  [3].  There  is  an  enormous 
amount  of  published  information  about  FDMs  in  various  scholarly  journals  and 
books.  This  section  describes  only  the  basics  of  FDMs,  and  the  interested  reader 
is  directed  to,  for  instance,  Gustafsson,  Oliger,  and  Kreiss  [4],  Finite  difference 
methods  focus  on  approximating  the  derivatives  of  the  solution  directly  at  a  set  of 
points  in  a  domain.  In  terms  of  the  calculus  of  finite  differences,  we  are  looking  to 
approximate  the  derivatives  by  linear  combinations  of  the  function  values  along  a 
grid  of  points  [5]. 

One  method  of  constructing  the  discretization  is  accomplished  by  a  Taylor  series 
expansion  on  a  selected  set  of  equally  spaced  points  (i.e.,. . .  <  x\-\  <  x\  <  x;+ \  <...). 
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For  instance,  suppose  you  wish  to  approximate  the  derivative  of  a  function  evalu¬ 
ated  at  Xi  using  grid  points  x;_ \  and  X{+\.  By  conducting  a  Taylor  series  expansion 
of  the  neighboring  values  around  X;,  we  can  construct  a  linear  combination  to  get 
an  accurate  approximation  of  the  derivative. 


For  example,  let  us  consider  the  one-dimensional  advection  equation. 


du  dll  _ 

~dt  +  dx  =  °' 


(1.1) 


where  u  =  u(x,t )  is  the  solution.  We  assume  an  appropriate  set  of  initial  conditions 
and  boundary  conditions  for  u(x,t).  Using  the  center  stencil,  we  can  Taylor  expand 
in  space  around  X{,  with  grid  spacing  Ax,  to  find  a  derivative  approximation  for 
where  U{  =  u(x,,t)  [3].  Starting  with 


OUi  AX  <7  Mf 

Mj+1  =  Ui  +  AX—  +  —  +  0(Ax  ),  (1.2) 

dll,  A  x2d2Ui  o, 

Hi- 1  =  Ui  -  AX—  +  —  +  0( AX6),  (1.3) 

(Uj+ 1  and  Ui- \  are  not  boundary  points)  and  taking  the  linear  combination  of  (1.2) 
and  (1.3)  yields 


dll{  Ui+i  Ui — i  ,  9 ,  ,  , 

*=^I^+0(a^'  (L4) 

which  can  be  substituted  into  (1.1)  to  yield  a  system  of  ordinary  differential  equa¬ 
tions.  Lastly,  we  can  solve  these  ODEs  computationally  through  a  numerical 
integration  scheme,  such  as  Runge-Kutta,  where 


dUi  _  dill  ^  Ui+1  lli- 1 
dt  dx  2AX 


Another  method  for  constructing  difference  formulas  is  to  build  a  polynomial 
interpolant  through  the  given  grid  of  points,  such  as  using  Lagrange  polynomials, 
and  evaluating  the  derivative  of  this  polynomial.  For  instance,  using  the  grid  of 
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points  [xi-i,Xi,Xj+i],  we  have  the  interpolant 

l 

Pi(x)  =  Yj  ui+kh(x),  (1-6) 

k=- 1 

where  L^{x)  is  the  Lagrange  polynomial  (see  Section  2.1).  The  derivative  approxi¬ 
mation  is  then 


duj 

dx 


1 

Lk(Xi)ui+k 

k=- 1 


ui+ 1  Uj- 1 

2AX 


(1.7) 


which  is  substituted  back  into  (1.1).  The  Taylor  series  approach  and  the  polynomial 
approach  are  equivalent  when  the  maximum  number  of  error  terms  are  eliminated 
from  the  linear  combination. 


1.2  Finite  Volume  Method 

Finite  volume  methods  differ  from  FDMs  in  that  instead  of  seeking  pointwise 
solution  values,  we  are  seeking  the  average  values  over  the  elements  using  approx¬ 
imations.  Depending  on  the  dimension  of  the  problem,  examples  of  elements  are 
cubes,  triangles,  or  intervals,  all  organized  in  an  unstructured  fashion  across  the 
physical  domain  [1],  From  the  average  values,  the  solution  is  reconstructed  in  or¬ 
der  to  evaluate  a  numerical  flux  to  tie  neighboring  elements  together  and  produce 
an  approximation  for  the  entire  system.  There  are  numerous  FVMs  and  here  we 
only  discuss  the  very  basics;  interested  readers  are  directed  to  LeVeque  [6]  for  more 
information. 

Let  us  again  consider  (1.1)  on  a  domain  that  is  represented  by  a  collection  of 
elements.  For  this  example,  let  us  use  elements  on  a  one-dimensional  domain  and 
each  element  is  denoted  by  the  index  i.  Let  us  further  define  a  different  set  of 
notation  with  ^  =  ut  and  jk  —  Ux_  Therefore,  (1.1)  becomes  lit  =  ~ux-  Generally,  the 
average  value  of  u.  on  an  element  is 

i  rR- 

Qi(t)=—  u(x,t)  dx  (1.8) 

Ax  Jl, 
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where  Ax,  in  this  case,  is  the  distance  between  the  left  and  right  boundaries  of 
the  one-dimensional  interval  and  Q/(f)  is  the  spatial  average  value  of  u  on  the  ith 
element.  If  we  were  using  cube-shaped  elements,  then  Ax  would  be  the  area  of  the 
cubic  element.  If  we  take  the  derivative  of  (1.8)  with  respect  to  time,  we  find 


dQi 

dt 


j_  rRi 

a*Jl; 


VLf  d-X 


(1.9) 


and  since  nt  =  —ux,  we  can  substitute  -ux  to  obtain 


dQi 

dt 


_j_  rRi 

A*  Jl, 


u  v  dx 


(1.10) 


and  integrate  (1.10)  to  find  an  equation  to  build  an  approximation  to  the  system: 


-jf  =  [u(Ri,t)  -  u(Li,t)] . 
dt  Ax 


(1.11) 


Equation  (1.11)  is  an  update  equation  for  the  average  value  of  the  solution  in 
the  ith  element,  where  the  time  derivative  of  the  average  value  changes  through 
the  fluxes  of  the  left  and  right  boundaries  of  the  element.  There  are  a  variety  of 
numerical  flux  methods  to  accomplish  this  approximation.  One  method  is  the  first 
order  upwinding  method,  meaning  that  we  select  the  average  value  of  the  interval 
always  on  the  "upwind"  side.  This  topic  is  discussed  in  further  detail  in  Chapter  3. 
Ultimately,  the  FVM  uses  the  average  element  values  from  across  the  domain  to 
reconstruct  an  approximation  for  the  system. 


1.3  Finite  Element  Methods 

Just  like  FDMs  and  FVMs,  there  are  many  different  FEMs  used  to  solve  PDEs. 
Similar  to  the  FVM,  the  physical  domain  is  mapped  into  a  reference  domain  of 
various  sized  geometrically  flexible  elements  known  as  the  grid.  After  the  grid  is 
built,  one  of  the  many  different  FEMs  is  applied  to  solve  the  system.  As  discussed 
at  the  beginning  of  this  chapter,  the  FEM  that  is  the  focus  of  this  paper  is  known  as 
a  Galerkin  method.  There  are  two  main  categories  of  Galerkin  methods,  Bubnov- 
Galerkin  and  Petrov-Galerkin  [7],  With  Petrov-Galerkin,  the  test  functions  and 
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basis  functions  are  different,  where  with  Bubnov-Galerkin,  the  test  functions  and 
basis  functions  are  the  same  and  reside  in  the  same  space  [1,  7],  We  discuss 
basis  functions  and  test  functions  further  in  later  chapters.  For  this  paper,  we  be 
focusing  on  Bubnov-Galerkin  methods,  often  called  simply  Galerkin  methods  [7], 
As  discussed  before,  the  two  main  categories  of  Galerkin  methods  are  Continuous 
Galerkin  (CG)  and  Discontinuous  Galerkin  (DG).  Both  methods  use  an  elemental 
approach  to  solving  the  system  in  integral  form,  but  the  difference  between  DG 
and  CG  is  mainly  whether  the  continuity  between  the  elements  is  enforced  strongly 
or  weakly  In  the  following  two  paragraphs,  let  us  lightly  touch  on  both  methods 
to  give  the  reader  a  small  insight  into  each  method  before  focusing  the  rest  of  the 
paper  entirely  on  element-based  DG. 

When  using  CG,  we  build  a  computational  domain  subdivided  into  various-sized 
elements,  similar  to  FVM.  However,  within  each  element  we  use  specially  selected 
degrees  of  freedom.  For  example,  in  the  one-dimensional  sense,  the  degrees  of 
freedom  could  be  a  grid  of  points  across  the  cell.  There  are  many  different  ways 
to  build  this  grid  of  points,  such  as  using  Legendre-Gauss  points,  equally  spaced 
points,  as  discussed  in  FDM,  or  some  other  combination  or  method.  For  the  purpose 
of  this  research,  we  use  Legendre-Gauss-Lobatto  (LGL)  points.  These  are  discussed 
in  further  detail  in  later  chapters,  but  they  are  a  set  of  points  that  cluster  toward 
the  boundaries  of  each  element  and  include  points  at  the  boundaries.  Discretizing 
the  integral  form  of  this  equation  yields  a  semi-discrete  scheme  with  elemental 
test  functions  [1].  These  test  functions  are  continuous  across  the  domain  in  CG 
since  we  enforce  the  element  boundary  points  to  be  equal  between  neighboring 
elements  [3].  Each  element  has  its  own  set  of  LGL  points;  however,  using  CG 
enforces  the  solution  at  the  element  boundaries  to  be  continuous  and  this  yields  a 
global  mass  matrix.  This  global  mass  matrix  is  for  the  entire  domain  and  must  be 
inverted  to  solve  time-dependent  problems;  this  can  be  computationally  expensive, 
depending  on  the  type  of  problem  [1]. 

Generally,  DG  is  similar  to  CG  in  requiring  numerical  interpolation  and  integration 
as  well  as  building  the  same  computational  domain;  however,  one  of  the  main 
differences  comes  from  the  elemental  boundaries.  In  CG,  we  require  the  boundaries 
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to  be  continuous,  while  in  DG  they  are  discontinuous  and  the  continuity  is  enforced 
weakly  The  domain  in  DG  is  still  represented  by  a  collection  of  elements;  the 
union  of  these  elements  is  accomplished  through  a  numerical  flux  similar  to  FVMs. 
Discontinuous  Galerkin  still  uses  the  same  space  of  basis  and  test  functions,  similar 
to  FEMs,  and  each  element  boundary  has  its  own  set  of  degrees  of  freedom  [1]. 
Therefore,  the  solution  is  typically  represented  by  a  set  of  piecewise  polynomials 
that  are  discontinuous  at  the  element  boundaries.  The  numerical  flux,  which  is 
discussed  in  greater  detail  later,  resolves  this  discontinuity  to  assist  in  finding 
the  final  solution.  Furthermore,  the  mass  matrix  is  constructed  locally  instead  of 
globally  and  this  allows  it  to  be  inverted  at  a  reduced  computational  cost,  yielding 
a  semidiscrete  scheme  that  is  explicit  [1].  Discontinuous  Galkerin  is  the  method 
used  for  the  following  study. 
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CHAPTER  2: 

Discontinuous  Galerkin  Foundation 


In  this  chapter,  the  information  and  terminology  maybe  confusing  to  the  reader 
who  has  no  knowledge  of  Galerkin  methods.  However,  just  like  pulling  out  a 
road  map  in  a  foreign  country,  the  intent  is  to  show  a  DG  road  map  and  apply 
these  concepts  to  a  one-dimensional  and  a  two-dimensional  problem  in  subsequent 
chapters.  Initially,  this  road  map  may  speak  a  different  language  to  the  reader,  as 
expected  for  anyone  new  to  this  numerical  method;  however,  by  the  end  of  this 
chapter  the  intent  is  for  the  DG  foundation  to  be  clear  and  understandable. 

2.1  Interpolation 

In  order  to  fully  implement  DG,  we  first  need  to  construct  the  building  blocks. 
This  leads  to  polynomial  interpolation  [3].  Nodal  interpolation  is  the  construction 
of  an  Nth  degree  polynomial,  f^{x),  that  is  equal  to  a  given  function,  f(x),  when 
evaluated  at  a  set  of  jq  points,  with  i  =  0 that  is  _/jv(jq)  =  f(X[)  [7].  There  are 
many  different  interpolation  methods  and  our  focus  for  this  section  is  Lagrange 
interpolation. 

Let  us  explain  Lagrange  interpolation  through  an  example.  We  are  going  to  ap¬ 
proximate  the  following  Gaussian  function  in  one  dimension, 

f(x)  =  e~Ax2r  x  G  [—1,  +1].  (2.1) 

We  begin  by  generating  a  grid  of  LGL  points  across  the  given  domain  in  one 
dimension.  As  discussed  in  Section  1.3,  LGL  points  are  a  specially  selected  set  of 
points  that  includ  the  boundaries  and  cluster  toward  the  ends  of  the  domain.  This 
clustering  of  points  toward  the  ends  of  the  domain  assists  in  avoiding  the  Runge's 
phenomenon  during  the  approximation,  which  is  the  oscillation  of  an  interpolation 
near  the  boundaries.  This  phenomenon  is  evident  when  equally  spaced  grid  points 
are  used  with  high  polynomial  orders.  Therefore,  clustering  points  toward  the  ends 
of  the  domain  helps  improve  the  interpolation  [8].  After  the  grid  of  LGL  points  is 
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generated,  we  construct  the  Lagrange  basis  functions,  L,(x),  by  the  equation 

t — r  (pC  ~  X ;) 

ll.v,  v.i-  i  =  0' . N'  {22) 

i*i 

where  X{  are  the  LGL  points  defined  on  the  domain  [ — 1, +1].  Using  (2.2)  for  all 
the  LGL  points  in  the  domain  generates  a  set  of  Nth  order  polynomial  basis  func¬ 
tions  associated  with  each  point  X[  in  the  domain.  With  these  basis  functions,  the 
approximation  /n(x)  is 


N 

=  (2.3) 

i= 0 

where /,  =  f(Xj)  for  i  =  0 N.  Moreover,  since  (2.2)  has  the  cardinal  property 


Li(xj)  — 


for  i  =  j 
for  i  ±  ]  ' 


(2.4) 


we  have  the  interpolation  property  /;  =  /v(-b)-  In  the  end,  f\j(x)  is  an  Nt]l  order 
polynomial  approximation  for  f(x). 


We  now  return  to  finding  an  approximation  for  representing  (2.1).  Figure  2.1 
shows  the  polynomial  interpolation  of  (2.1)  using  N  =  2,  4,  8,  and  16  with  N  +  1 
LGL  interpolation  points.  In  Figure  2.2,  we  consider  the  error  norms 


e  Wli 


e  11/2 


z 

£iI/n(4)-M)I 

1  /(4)  1 

\ 

X£i(/n(4)-/(4))2 

LfUWh))2 

maxl<k<50  1  /n(4)  ~f(Xk)  1 

max^^o  |  f(xk)  | 

(2.5) 

(2.6) 

(2.7) 


where  we  have  used  an  equally  spaced  grid  of  50  points,  x^,  within  the  domain 
[— 1,+1].  Figure  2.2  displays  the  error  norm  of  N  =  1  through  N  -  40  polynomial 
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interpolations. 


As  you  can  see  in  Figure  2.1,  the  16th-  order  polynomial  interpolation  does  a  good 
job  approximating  the  function  as  compared  to  the  exact  solution  (denoted  by  *  to 
differentiate  from  the  various  interpolations).  In  Figure  2.2,  we  show  the  error  in 
the  interpolation  vs.  the  Nth  order  polynomial  interpolant  based  off  of  (2.5),  (2.6), 
and  (2.7).  As  the  figures  show,  the  approximation  gets  to  machine  precision  around 
the  32nd  order  polynomial.  Lastly,  the  reason  that  Figure  2.2  has  various  plateaus 
within  the  convergence  rate  is  because  f(x)  is  an  even  function. 


Figure  2.1:  Lagrange  Interpolation  ofe~ixl  at 
different  sets  of  LGL  points  for  various  poly¬ 
nomial  orders  ofN. 


Figure  2.2:  Approximation  Error  based 
off  (2.5),  (2.6),  and  (2.7)  measured  at  poly¬ 
nomial  orders  N  =  1  through  N  =  40. 


2.2  Integration 

The  next  building  block  for  DG  is  the  ability  to  conduct  numerical  integration. 
When  performing  numerical  integration,  also  known  as  quadrature,  the  analysis  is 
similar  to  that  of  interpolation.  Integrating  the  Lagrange  interpolation  (2.3)  gives 
the  integral  approximation 

X+i  r*+i  ^  r'+i 

fN(x)dx=  YLi(x)fidx  =  Yf(xi)  Li(x)dx.  (2.8) 

i  J-i  7T  J-i 
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Defining  the  quadrature  weights 


f+1 

W(  =  J  Li(x)dx 

gives  the  numerical  integration  method  [8] 


(2.9) 


p+l  N 

I  fN(x)dx  =  ^  Wifi. 

^-1  z=0 


(2.10) 


These  quadrature  weights  only  have  to  be  precomputed  once  with  the  LGL  points 
and  not  at  each  evaluation  of  a  particular  integral  [8].  Figure  2.3  is  the  comparison 


Figure  2.3:  Numerical  Integration  Error  evaluated  with  the  Euclidean  Norm  vs.  the  Nth  order 
quadrature. 

of  the  numerical  integration  approximation  to  the  actual  value  of  1  4x^dx  ~ 
0.882081.  The  Euclidean  Norm  is  used  for  the  calculation  of  the  error  due  to  the 
integration  of  f(x)  producing  a  scalar.  As  you  can  see,  the  error  decreases  to 
machine  precision  around  the  18^'  order  quadrature. 
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In  comparing  Figure  2.2  to  Figure  2.3,  we  see  that  when  using  the  same  order 
polynomial,  integration  is  more  accurate  than  interpolation.  In  general,  using 
N  + 1  points,  we  can  construct  an  Nth  order  polynomial.  Since  we  fixed  the  LGL 
points  at  the  boundaries,  -1  and  +1,  there  are  only  N—  1  LGL  points  left  to  choose. 
We  also  have  N  + 1  quadrature  weights  to  choose.  Therefore,  we  have  N  —  l  points 
and  N  + 1  weights,  which  means  there  are  2 N  degrees  of  freedom.  We  can  thus  set 
the  degrees  of  freedom  so  that  we  integrate  a  polynomial  of  order  2 N  - 1  exactly 
(as  this  polynomial  has  2 N  coefficients)  [3].  For  example,  we  can  exactly  integrate 
a  fifth  degree  polynomial  using  N  + 1  =  4  LGL  points;  even  though  we  can  only 
construct  a  third-order  interpolating  polynomial  using  these  points. 


2.3  Concept  of  the  Mass  Matrix 

After  discussing  interpolation  and  integration,  let  us  now  consider  the  concept 
of  the  mass  matrix.  The  mass  matrix,  denoted  by  M,  is  used  for  integration. 
We  discuss  the  integration  of  the  product  of  two  interpolants  exactly  and  inexactly. 
Using  exact  integration  results  in  a  full  mass  matrix,  while  using  inexact  integration 
results  in  a  diagonal  mass  matrix.  For  example,  suppose  we  wanted  to  integrate 
f(x)  with  a  test  function,  i p(x), 

J'  f(x)ip(x)dx,  (2.11) 

over  the  domain  of  a  single  element  from  [-1,+1].  Test  functions  are  defined  in 
further  detail  in  Section  3.3,  but  for  this  section  all  we  need  to  consider  is  that  ip(x) 
is  a  function  that  resides  in  the  same  space  as  f(x). 

Suppose  both  i p(x)  and  f(x)  can  be  numerically  approximated  with  Lagrange  poly¬ 
nomials.  Then 


N 

ip(x)  ~  1 Pn(x)  =  Y  xf>(xi)Li(x)  and  (2.12) 

i= 0 
N 

f(x)  «  /n(x)  =  Yj  fixi)Li(x).  (2.13) 

i= 0 


11 


This  leads  to  the  integral  approximation 


X+i  /^+i 

f(x)ip(x)dx  x  I  f^(x)ip^(x)dx  =  ££/(*i)[L  i(x)Lj(x)]ip(xj)dx, 

1  ^_1  ;=0  ;'= 0 


which  is  further  simplified  to 

N  N 


N  N  r  ^+1 

EI>W>  J .  1 

V— n  v— n 


Li(x)L;(x)dx 


i= 0  j=0 


/(x,)  =  xpTMef, 


(2.14) 


where 


Li(x)Lj(x)dx  =  M-y 


(2.15) 


and  Me-j  is  a  full  matrix  per  element  that  can  integrate  an  order  2 N  function  exactly. 
Furthermore,  ij>  and  /  are  the  evaluation  of  ip(x)  and  f(x)  at  the  LGL  points  stored 
in  column  vectors: 


'/(*  o)' 

\p(x0) 

fix  l) 

/  4  = 

Vix  l) 

/(*n). 

.i/^n). 

Inexact  integration  is  similar  to  exact  integration.  Instead  of  numerically  approxi¬ 
mating  the  functions  separately  with  Lagrange  polynomials  as  seen  above,  we  first 
multiply  the  function  and  test  function  together  at  the  LGL  points  and  produce  a 
new  interpolant  for  numerical  integration: 


I  fN(x)ipN(x)dx  ~  I  Yf{xiMxi)Li{x)dx=YMm.  (2.16) 

to  to 


Once  again  L/(x),  under  integration,  generates  the  weights,  w-u  as  seen  in  (2.9), 
and  Wi  is  only  computed  once  at  the  LGL  points.  These  weights  can  be  stored 
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in  a  diagonal  mass  matrix,  called  M,  and  used  in  an  one-dimensional  problem. 
Therefore,  (2.11)  can  now  be  written  in  the  matrix  form  below: 

r+1 

J  f(x)xp(x)dx  ~  i[)TMf.  (2.17) 

As  discussed  in  Section  2.2,  using  this  method  with  LGL  points  allows  us  to  inte¬ 
grate  a  polynomial  order  of  2 N  - 1  exactly.  Throughout  the  remainder  of  this  thesis, 
we  only  consider  inexact  integration.  This  is  because  the  diagonal  structure  of  the 
mass  matrix  greatly  simplifies  the  implementation  of  the  method.  Additionally,  the 
Jacobians  and  surface  Jacobians  are  easier  to  handle  when  moving  into  multiple 
dimensions.  We  discuss  more  about  Jacobians  and  surface  Jacobians  in  subsequent 
sections  and  chapters. 


2.4  Concept  of  the  Differentiation  Matrix 

The  differentiation  matrix  has  a  purpose  similar  to  that  of  the  mass  matrix,  but  it  is 
obviously  used  for  differentiation.  Like  the  mass  matrix,  the  differentiation  matrix 
is  found  using  Lagrange  polynomials  by 


dLj(xi) 

D'i  dx 

dfN(Xi)  _  dLj(xj )  r  ^  r 

dx  ~  dx 

i= 0  i= 0 


(2.18) 

(2.19) 


where  it  should  be  evident  from  (2.19)  that  Df  approximates  ^  at  all  the  nodes. 
As  you  can  see,  D  is  simply  the  derivative  of  the  Lagrange  polynomials  evaluated 
at  the  LGL  nodes.  You  should  further  notice  that  D  is  a  full  matrix. 


Let  us  examine  how  the  differentiation  matrix  is  used  for  discretization.  Consider 
the  integral 


df(x) 

dx 


i p(x)dx, 


(2.20) 


where  xp(x)  is  a  test  function  that  resides  in  the  same  space  as  f(x).  Once  again,  the 


13 


functions  are  numerically  approximated  with  a  set  of  Lagrange  polynomials 


N 

$n(x)  ~  ^(Xi)Li(x), 

1=0 

dfN(x)  a  \dU(x) 

i=  0 


(2.21) 

(2.22) 


leading  to 


df(x ) 
dx 


ip(x)dx 


dLi(x) 

dx 


i p(xj)dx, 


which  further  simplifies  to 


(2.23) 


dLi(x) 

dx 


L;(x)dx 


/(x^/MD/. 


(2.24) 


As  you  can  see,  (2.18)  is  nestled  right  in  the  middle  of  (2.24).  Once  again,  both  i/> 
and  /  are  the  evaluation  of  ip(x)  and  /(x)  at  the  LGL  points. 


2.5  Change  of  Variables 

Mapping  from  a  physical  space  to  a  computational  reference  space  requires  a 
change  of  variables.  In  this  case,  we  are  going  to  map  from  x  e  [xn,x„+i]  to  £  e 
[— 1,+1],  where  n  =  0 ,...,N.  Furthermore,  a  Jacobian  arises  when  conducting  this 
change  of  variables  [9].  In  one  dimension,  the  Jacobian  is  simply  |,  where  the 
element  size  is  h  =  Ax  for  an  equally  spaced  grid  of  elements.  This  one-dimensional 
Jacobian  is  used  extensively  in  Chapter  3  for  discretization  and  is  annotated  by 
/  =  |.  In  two  dimensions,  it  is  a  little  more  complicated  and  is  discussed  in  greater 
detail  in  Chapter  4.  However,  the  concept  remains  the  same  for  one  dimension  and 
two  dimensions. 

In  one  dimension,  we  want  the  computational  reference  domain  to  be  from  £  e 
[-1,  +1]  for  two  points  from  a  linear  element,  xn  and  xn+\.  This  concept  is  depicted 
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in  Figure  2.4. 


Xn  Xn+i  =  _1  g1  =  +1 

Figure  2.4:  Physical  Space  to  Computational  Space 


The  change  of  variables  for  Figure  2.4  is  obtained  by  a  linear  combination 


*(£)  = 


(1+5)' 


Xfl+lr 


(2.25) 


that  allows  us  to  approximate  the  coordinates  of  the  element.  For  example,  if  E,  -  —1 
in  (2.25),  then  that  maps  to  xn  and  if  t  =  +1,  that  maps  to  xn+\ .  Taking  the  derivative 
of  (2.25)  yields 


dx  _  (xn+i  -xn)  _  h 

H =  2  =  2' 


(2.26) 


and  dx  =  | dE,  [3].  Therefore,  using  (2.25)  and  (2.26)  for  the  change  of  variables  for 
integration  gives 


rXn+ 1  r+1  h 

J  f(x)dx  =  J  i  -f(x(£))dt,  (2.27) 

where  xn  and  xn+i  are  the  left  and  right  boundaries  for  the  physical  element. 
Constructing  this  change  of  variables  for  each  element  is  a  key  foundation  for  local 
element  based  Galerkin  methods.  Instead  of  building  different  matrices  for  each 
element  and  solving  them  individually,  we  can  now  construct  the  required  matrices 
for  a  reference  element  and  then  use  the  metrics  terms  to  scale  the  reference  element 
to  the  physical  space  [3].  Metric  terms  are  discussed  in  greater  detail  in  Chapter 
4.  As  for  using  the  differentiation  matrices,  D,  the  change  of  variables  requires  the 
chain  rule,  which  creates  a  =  ft  multiplied  by  ||  =  cancelling  these  terms  out. 
The  equation  below  is  a  visual  example  for  the  above  sentence  showing  the  change 
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of  variables  for  (2.20): 


m^dx- 

dx 


’dx  dx 


df{xm 


tpTMDf . 


This  section  assists  in  the  construction  of  the  one-dimensional  grid  further  investi¬ 
gated  in  Section  3.1. 
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CHAPTER  3: 

One-Dimensional  Discontinuous  Galerkin 


Let  us  start  working  through  the  Discontinuous  Galerkin  method  for  the  linear 
acoustic  wave  equation  in  one  dimension.  Consider  the  equation 

=  c2V2p,  xeQc  IRd,  t  >  0,  (3.1) 

where  p  is  the  pressure,  c  >  0  is  the  wave  speed,  and  d  is  the  spatial  dimension.  On 
the  boundary,  represented  by  dD,,  we  impose  periodic  boundary  conditions.  We 
are  now  going  to  split  this  second-order  wave  equation  into  two  equations  with 
the  auxiliary  variable,  co,  equalling  the  time  derivative  for  pressure  [2]: 

<3-2> 

!="•  (3-3) 

3.1  One-Dimensional  Grid 

Figure  3.1  shows  an  example  of  a  one-dimensional  grid  of  equally  spaced  elements, 
with  polynomial  order  N  =  6  and  three  elements,  covering  Q  =  [— 1,+1],  The  LGL 
points,  depicted  by  red  x  symbols,  are  used  for  the  interpolation  points  per  element. 
Polynomial  order,  N,  is  defined  as  the  maximum  order  of  a  polynomial  that  can  be 
represented  exactly  on  each  element.  Obviously,  the  grid  varies  depending  on  N 
and  number  of  elements  (Ne).  As  you  can  see,  the  LGL  points  are  not  evenly  spaced 
across  the  individual  elements  and  cluster  towards  the  boundaries  of  each  element. 
There  are  N  + 1  LGL  points  per  element,  where  each  element  is  represented  by  Q;, 
with  j  =  1  ,...,Ne.  Remember,  this  is  a  one-dimensional  grid  with  Q  =  [— 1,  +1] 
and  each  element  (Q y)  is  mapped  to  E  e  [ — 1, +1]  through  a  change  of  variables,  as 
discussed  in  Section  2.5.  This  one-dimensional  grid  is  simple;  however,  in  two 
dimensions,  the  grid  is  a  little  more  complicated  and  is  discussed  in  Chapter  4. 
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Figure  3.1:  One-Dimensional  Grid  with  N  =  6  and  Ne  =  3. 


3.2  Variational  Form 

Understanding  (3.2)  and  (3.3)  motivates  the  variational  form  that  is  needed  for  us 
to  build  the  set  of  equations  for  DG.  This  section  explains  the  process  of  finding 
the  variational  form  of  (3.2)  and  (3.3).  Furthermore,  the  discretization  of  these 
equations  yields  the  necessary  equations  for  the  DG  numerical  approximation. 
First,  we  need  to  focus  on  finding  a  variational  form. 

3.2.1  First  Variational  Equation 

Equation  (3.4)  is  the  most  trivial  of  the  three  variational  form  equations.  If  ^  =  co, 
it  follows  that 


Using  (3.4)  in  Section  3.4  assists  in  finding  the  needed  equations  for  the  DG  ap¬ 
proximation. 


3.2.2  Second  Variational  Equation 

Given  that  ^  =  c2V2p,  by  multiplying  this  equation  by  a  test  function,  1 p,  and 
integrating  over  Q y  yields 


We  define  test  functions  in  further  detail  in  the  following  section.  By  conducting 
integration  by  parts  on  (3.5)  and  inserting  a  numerical  flux  at  the  boundary  of  the 
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element,  we  are  able  to  find  the  second  equation  in  the  variational  form  to  be 


+  c2Vip-Vp 


f 

J  d  Q, 


ip  n  ■  (Vp)* 


(3-6) 


For  notation  purposes,  a  *  term  in  the  equation  identifies  the  numerical  flux.  Using  a 
numerical  flux  is  a  numerical  technique  for  coupling  the  elements  that  is  commonly 
used  in  DG.  For  example,  ( Vp )*  is  the  coupled  numerical  flux  computation  for  the 
gradient  of  the  pressure  at  the  boundaries  (dQ, j)  for  neighboring  elements  within 
the  domain.  An  explicit  computation  of  the  numerical  flux  is  discussed  in  greater 
detail  in  Section  3.5. 


3.2.3  Third  Variational  Equation 

In  order  to  find  the  last  variational  equation,  we  are  going  to  multiply  the  gradient 
of  f  ^  -  co \  with  the  gradient  of  a  test  function  (Vxp): 


f  Vi/>V 
Jo, 


(3.7) 


Let  us  focus  on  the  last  portion  of  (3.7).  By  conducting  integration  by  parts  we  find 


r  r 

VxpVco  = 

JQ; 


(Vi p-n)co- 


(3.8) 


By  introducing  a  numerical  flux  to  the  boundary  term  of  (3.8)  and  substituting  it 
into  (3.7),  we  find  that 


(Vi p-n)co*+  f  (vV) 
JdQj  JQ;  ; 


co  =  0. 


(3.9) 


Let  us  focus  now  on  the  last  term  from  (3.9).  Conducting  integration  by  parts  again 
on  the  last  portion  of  (3.9)  yields 


f  (vy)o.= 

JQ;  J<9Q; 


(Vi p-n)co- 


n 

Vi/^Vcu. 

JQ; 


(3.10) 
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Next,  we  substitute  (3.10)  back  into  (3.9)  and  combine  like  terms  to  conclude  the 
following: 


f  Vi pV 

jQj 


dp 

Tt 


WxpVco  =  0. 


Finally,  from  (3.11),  we  arrive  at  the  final  variational  form 


(3.11) 


f  Vi//V 

JQj 


1  dp 
i  dt 


“")  =  f 

1  Jd  Qj 


(Vip  -n)(a>*  -co) . 


(3.12) 


Overall,  equations  (3.4),  (3.6),  and  (3.12)  together  are  the  variational  form  of  (3.1) 
that  we  discretize  using  DG. 


3.3  Test  Functions 

Using  Discontinuous  Galerkin,  we  are  looking  for  p  and  co  from  (3.2)  and  (3.3) 
that  satisfies  the  variational  form  found  in  Section  3.2  for  all  piecewise  polynomial 
test  functions.  Since  we  are  working  to  build  an  approximation  to  a  function  with 
piecewise  polynomials,  our  space  is  the  space  of  Nth -degree  piecewise  polynomials, 
with  the  objective  of  solving  for  a  piecewise  polynomial  that  represents  the  solution. 
For  example,  suppose  we  wanted  to  find  a  constant  approximation  for  f(x),  such 
that  the  error  was  orthogonal  (in  an  integral  sense)  to  all  other  constants.  That  said, 
we  want  to  find  a  G  R  such  that 


X+l  p+L 

ip(x)f(x)dx  =>  J  (C)(a  -  f(x))dx  =  0, 
for  all  C  G  1R.  If  f(x)  =  x2,  then  the  problem  is  to  find  a  such  that 

J+1,_  1  1+1 


(Ca  -  Cx  )dx  = 


Cax--Cx 


=  0, 


which  is  simplified  to 


2Cn -  =  2 c(a-^-)  =  0. 


(3.13) 
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Therefore,  a  =  ^  is  an  approximation  for  x2  found  using  a  0th  degree  test  function 
(C).  Expanding  this  concept  to  DG  is  essentially  the  same,  except  we  are  using  Nth- 
degree  piecewise  polynomials  as  the  test  functions  to  find  a  Nth  degree  polynomial 
approximation  for  the  solution  on  each  element.  We  use  this  concept  for  the 
remainder  of  this  paper  and  Section  3.3.1  further  highlights  the  use  of  test  functions 
for  the  discretization  in  Section  3.4  and  Section  4.2. 

3.3.1  Test  Functions  in  Discretization 

As  discussed  in  the  previous  section,  we  are  using  Nth- degree  piecewise  polynomi¬ 
als  as  test  functions  for  DG.  In  the  following  section  and  chapter,  we  discretize  the 
variational  form  to  find  a  discrete  approximation  for  the  continuous  equation  in 
one  dimension  and  two  dimensions.  These  discrete  approximations,  with  respect 
to  ip,  all  have  the  similar  form  of 

ipT  (Ap  -  Bco)  =  0,  (3.14) 

where  A  and  B  are  matrices  and  p  and  co  are  the  solution  vectors.  Since  (3.14)  holds 
for  all  ip,  then  Ap  -  Bco  must  equal  zero.  This  concept  holds  true  for  Section  3.4 
and  Section  4.2.  Therefore,  in  what  follows,  ip  is  not  listed  in  the  final  discrete 
approximations. 

3.4  One-Dimensional  Discretization 

In  this  section,  we  derive  discrete  versions  of  Equations  (3.4),  (3.6),  and  (3.12).  The 
idea  of  discretization  is  to  replace  a  continuous  equation  with  a  consistent  discrete 
approximation.  As  you  have  probably  noticed,  with  Galerkin  methods  we  use  the 
variational  form  for  discretizing  equations  [3].  We  focus  on  developing  a  spatial 
discretization  of  the  equations  and  time  is  integrated  separately  with  a  Runge-Kutta 
method.  We  isolate  the  time  derivatives  and  combine  the  discrete  forms  of  (3.4) 
and  (3.12)  into  one  discrete  equation  [10].  As  discussed  in  Section  2.5,  the  change 
of  variables  requires  the  Jacobian  of  /  =  |  for  an  equally  spaced  grid  of  elements. 
As  a  reminder,  the  concept  from  Sections  3.3.1  is  applied  to  the  discretization  in 
Section  3.4.2  and  Section  3.4.3.  We  discretize  each  equation  individually. 
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3.4.1  First  One-Dimensional  Discretization  Equation 

Using  the  concepts  established  in  Chapter  2,  the  discrete  version  of  (3.4)  is 

lTMJ^-lTMJa)  =  0.  (3.15) 

at 

Within  (3.15),  1  is  a  vector  of  ones  and  is  needed  because  in  (3.4)  we  are  integrating 
against  the  function  xp(x)  =  1  [10].  Additionally,  ^  and  co  are  column-vector  ap¬ 
proximations  of  the  solution  at  the  grid  points.  This  is  similar  to  /  from  Section  2.3. 


3.4.2  Second  One-Dimensional  Discretization  Equation 

The  following  is  the  discrete  version  of  (3.6): 


M/^  +  c2rrDTMDp  =  c2rl 


'  r)p \ 


'  r)p ' 


ot 


(3.16) 


The  column  vectors  e n  and  cq  consist  entirely  of  zeros  except  for  the  last  "row" 
in  ej v  and  the  first  "row"  of  eo  being  ones.  Using  and  eo  is  a  way  of  ensuring 
that  the  first  and  last  portion  of  information  from  the  numerical  flux  of  |||  j  and 

||  j  is  used  for  the  calculation.  The  numerical  flux  is  discussed  in  greater  detail 
in  Section  3.5,  but  it  is  essentially  a  method  for  coupling  the  solution  on  either  side 
of  an  element  boundary.  Additionally,  discretization  of 


from  (3.6)  requires  two  differentiation  matrices;  DT  takes  the  derivative  of  the  test 
function,  M  is  the  integral,  and  D  is  the  derivative  of  p.  This  accounts  for  the  DTMD 
portion  of  (3.16). 


3.4.3  Third  One-Dimensional  Discretization  Equation 

The  following  is  the  discrete  version  of  (3.12): 


J~1DTMDd^  -  J~1DtMDo)  =  r1DTeN(cv*N-elja))-r1DTeo(a)*0-elcv).  (3.17) 
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Once  again,  the  Vi/>-v||  and  fQ  Vip-Vco  from  (3.12)  are  accounted  for  by  the 
DtMD  portions  and  co*  is  computed  through  the  flux  at  the  element  boundaries. 


3.4.4  Combination  of  One-Dimensional  Discretization  Equations 

Now  that  we  have  three  discrete  equations,  our  goal  is  to  get  the  problem  into  the 
form  of  two  equations  and  two  vector  unknowns.  Combining  (3.15)  and  (3.17) 
yields 

(/_1DrMD  +  1MJ)  ^  =  (j rlDTMD  +  1  Mj)  w  +  J  rDTeN  (a >*N  -  eTNco)  ^ 

-r1DTe0(co*0-eT(o). 


For  notational  purposes,  1  is  a  necessary  square  matrix  of  ones  of  size  (N  +  1,N  + 1) 
because  (3.15)  is  a  scalar.  Additionally,  letting  {j~1DTMD  +  1  Mj)  be  equal  to  the 
variable  Mi  and  simplifying  (3.18)  produces 

Mi  ^  =  Mtco  +  J~lDT (eNa)*N  -  e0a)*0)  -  /_1DT ( eNe^co  -  e0 ejm) .  (3.19) 

For  the  second  equation,  by  rearranging  (3.16)  we  can  isolate  ^  as  seen  in  the 
following: 


MJ 


dco 

dt 


-c2J~1DTMDp  +  c2/-1 


2n\ 


'dp' 

M, 


-Col 


N 


'dpS 


(3.20) 


Equations  (3.19)  and  (3.20)  are  the  two  systems  of  ordinary  differential  equations 
that  we  solve  numerically  using  a  Runge-Kutta  method. 


3.5  Numerical  Flux 

When  using  a  DG  method,  we  have  to  account  for  the  discontinuity  that  exists 
between  the  elements.  Meaning,  when  each  element  is  represented  by  a  polynomial 
for  co  and  p,  sampling  co  and  p  at  the  boundary  of  neighboring  elements,  as  shown  in 
Figure  3.2,  yields  different  results  [3].  The  numerical  flux  is  a  method  for  enforcing 
(approximately)  continuity  between  the  elements.  There  are  multiple  flux  methods, 
but  in  what  follows  we  only  discuss  the  central  and  upwinding  fluxes. 
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Figure  3.2:  One-Dimensional  Element  Boundary 


The  central  flux  is  the  easiest  method  to  understand,  but  it  does  not  necessarily 
produce  the  best  results.  The  central  flux  is  the  average  of  the  values  at  each 
boundary  element.  In  Figure  3.2,  the  neighboring  elements,  Q;  and  Clj+i,  have 
different  approximations  at  the  boundary  point  (BP).  For  notation  purposes,  the 
superscript  (+)  represents  the  left  element's  gridpoint  and  the  superscript  (-)  repre¬ 
sents  the  right  element's  gridpoint.  Furthermore,  as  in  Section  3.4,  the  superscript 
(*)  denotes  the  numerical  flux  term.  This  notation  comes  into  further  use  in  the 
following  section  and  chapters.  The  central  flux,  for  both  co  and  is  defined  by 
the  following: 


eo*=-(co  +co+), 
dp*  _  1  (dpY  +(fy\ 

dx  2  \dx)  \dx) 


(3.21) 

(3.22) 


Using  the  central  flux  is  a  useful  beginning  step  when  coding  the  flux  into  any 
algorithm.  The  central  flux  is  simple  and  still  achieves  convergence  when  coded 
correctly.  However,  central  fluxes  often  have  suboptimal  convergence  rates  and 
can  lead  to  oscillatory  approximations  due  to  a  lack  of  dissipation.  To  improve  this, 
we  need  a  numerical  flux  that  does  not  just  average  the  information  at  the  element 
boundaries,  but  can  account  for  the  physical  propagation  of  information  across  the 
boundaries.  This  leads  us  to  the  upwinding  flux. 

An  upwinding  flux  is  a  method  for  computing  a  flux  across  the  element  boundaries 
based  on  the  physical  propagation  of  information.  In  the  acoustic  wave  equation, 
the  waves  that  are  propagating  through  the  physical  system  can  move  in  both 
directions  across  the  element  boundaries.  We  are  looking  at  the  information  that 
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is  "upwind"  of  the  wave's  direction  of  motion,  hence  called  an  upwind  flux.  For 
example,  for  the  advection  equation,  if  the  wave  is  moving  across  the  element 
boundary  from  the  left  to  the  right,  we  select  the  information  from  the  left  element 
and  vice  versa  for  the  wave  moving  from  right  to  left  across  the  boundary.  However, 
for  the  wave  equation  we  want  to  decouple  the  wave  propagation  into  two  one-way 
advection  equations  to  take  advantage  of  this  "upwind"  concept. 

Since  we  are  working  in  one  dimension,  let  us  take  (3.2)  and  derive  the  upwind 
flux.  Defining  the  auxiliary  variable  q  =  px,  we  rewrite  (3.2)  as 

cot  =  c2qx.  (3.23) 
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We  define  a  new  set  of  variables  r\  and  known  as  the  characteristic  variables,  by 


r=V~1U  = 


n 

ri 


and  let  A  be  the  diagonal  matrix  of  eigenvalues 


[v_1AV]  =  A  = 


Ai 

0 


0 

A2 


We  now  have  a  system  of  first  order  one-way  advection  equations  that  allow  us  to 
upwind  because  the  propagation  of  the  waves  are  decoupled  across  the  element 
boundaries  by  the  characteristic  variables,  r\  and  ri,  as  seen  with 


rt  -  A  rx  =  0  => 


n 

Ai 

0" 

D 

'o' 

A 

t 

0 

A2 

r2 

X 

0 

(3.25) 


By  analysis  of  the  function  within  the  element,  with  f\  and  being  the  initial 
conditions  for  r\  and  r 2,  we  find  that 


r\(x,t)  =  f\(x  +  A\t)  =  fi(x  +  ct),  and  r2(x,t)  =  fi{x  +  A2O  =  fiix-ct), 

which  further  describes  the  propagation  of  the  wave  across  the  element  boundaries 
with  r\,  flowing  right  to  left,  and  rj_,  flowing  left  to  right.  Figure  3.3  depicts  this 
propagation.  We  know  the  direction  of  propagation  and  can  now  find  the  flux 
values.  By  letting  r*  =  r~  and  r*  =  ,  we  have  chosen  the  numerical  fluxes  as  the 


Q, 


R  Q 


7+1 


Figure  3.3:  One-Dimensional  Element  Boundary  Decoupled  Wave  Propagation 
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upwind  values.  Knowing  that  r  =  V  1U,  we  can  easily  find  a  system  of  equations 
to  solve.  Note  that 


n 

■  i 

2c 

r 

2 

CO 

— 

1 

1 

/2 

2c 

2 

T 

so  with  r*  =  r1  and  r*  =  we  have 


=  2cCV~  +  2CI~  = 


1  1 

— co*  +  -q*  and 

2c  2 7 

(3.26) 

1  1 

*  .  * 

—  ~~co  +  —q  . 

2c  2 7 

(3.27) 

We  now  have  two  equations  and  two  unknowns  where  we  can  solve  for  co*  and  q*: 


!)•  W(0>--O+i(i|++<r);  (3.28) 

">'  =  j  (ar  +  a>+)  +  ^  — 1/+) .  (3.29) 

Lastly,  with  q  =  (3.28)  and  (3.29)  are  the  upwinding  flux  equations  for  cu  and 

in  one  dimension.  It  should  be  further  noted  that  both  central  flux  equations,  (3.21) 
and  (3.22),  are  included  in  (3.28)  and  (3.29).  The  portion  of  (3.28)  and  (3.29)  that 
include  the  variable  c  is  considered  the  upwinding  portion  and  is  used  separately 
for  the  upwind  flux  energy  analysis  in  Chapter  5. 


3.6  One-Dimensional  Discontinuous  Galerkin  Results 

After  laying  the  groundwork,  we  now  apply  these  concepts  to  an  actual  problem 
in  one  dimension  using  an  upwind  flux.  We  use  the  exact  solution 

p(x,  t )  =  sin(n7ix)  sin(nnt ),  (3.30) 

dp 

—  =  co  =  nnsm(nnx)cos(nnt),  (3.31) 

dt 

on  the  domain  Q  =  [-1,  +1]  with  periodic  boundary  conditions  where  n  is  an  integer. 
To  build  the  initial  condition,  (3.30)  and  (3.31)  are  evaluated  at  t  -  0  across  a  one¬ 
dimensional  grid  of  LGL  points  where  x  £  [— 1,  +1].  These  initial  conditions  are  used 


27 


within  (3.19)  and  (3.20)  and  evaluated  through  a  Runge-Kutta  five-stage  fourth- 
order  (RK54)  accurate  iterative  method  [11].  Runge-Kutta  methods  are  numerical 
methods  that  approximate  a  system  of  ordinary  differential  equations.  With  the 
RK54  method,  (3.19)  and  (3.20)  are  evaluated  through  five  stages  per  time-step 
along  the  grid  of  LGL  points.  All  of  these  evaluations  are  combined  to  produce  a 
fourth  order  accurate  or  higher  approximation  [8].  The  following  section  discusses 
the  results  of  applying  a  DG  method  for  the  above  one-dimensional  problem.  The 
DG  implementation  can  be  found  in  Appendix  B. 

The  implementation  is  tested  using  various  polynomial  orders  and  numbers  of 
equally  spaced  elements  within  the  domain  listed  in  Table  3.1.  Once  again,  the 

Table  3.1:  Discontinuous  Galerkin  Tested  Information 


Polynomial  Orders  (N) 

Number  of  Elements  (Ne) 

N  =  4,6,8 

Ne  =  2,4,8,16 

N=  16 

Ne  =  2,4,8 

DG  algorithm  uses  inexact  integration  for  computational  speed  to  ensure  an  easily 
invertible  diagonal  mass  matrix.  Figure  3.4  shows  the  log  of  the  error  vs.  log  of  the 
number  of  elements  for  N  =  4,6,8  calculated  using  the  global  L2  error  norms.  As 
expected  for  both  co  and  p,  Figure  3.4  displays  increasing  convergence  rates  as  the 
polynomial  order  increases.  Increasing  the  order  of  the  local  approximation  gives 
the  fastest  convergence  rates,  as  apposed  to  increasing  the  number  of  elements, 
due  to  the  fact  that  the  global  error  is  dependent  on  the  polynomial  order  [1]  as 
depicted  by 


||  £  ||2  <  Chq.  (3.32) 

In  (3.32),  C  is  a  constant  that  is  not  dependent  on  the  element  size  h  but  does  depend 
on  the  final  time,  t,  of  the  solution.  In  this  case,  N  is  proportional  to  q  in  (3.32).  As 
the  polynomial  order  (N)  increases,  the  convergence  rates  increase. 

Table  3.2  displays  the  convergence  rates  for  the  tested  information  from  Table  3.1 
for  co  and  p.  As  you  can  see,  as  the  polynomial  order  increases,  the  convergence 
rates  increase.  Furthermore,  Table  3.2  depicts  the  convergence  rates  being  near  or 
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Figure  3.4:  Convergence  Rates  For  N  =  4,6,8  at  Ne  =  2,4,8, 16 


better  than  its  associated  polynomial  order,  where  the  convergence  rate  for  co  is  of 
order  N  and  the  convergence  rate  for  p  is  of  order  N  +  2. 

Most  DG  methods  are  expected  to  be  order  N,  N  +  \,  orN  +  1  [1].  However,  the 
results  yielded  N  and  N  +  2  when  using  the  same  space  of  functions  for  co  and  p. 
In  Appelo  and  Hagstrom's  paper  [2]  they  proved  that  when  p  is  an  N  + 1  order 
polynomial  and  co  is  an  N  order  polynomial  you  get  optimal  convergence  of  N  + 1 
for  this  method  [2].  In  Table  3.2  and  3.3,  these  convergence  rates  are  an  observation; 
a  proof  would  require  further  analysis. 
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Table  3.2:  Convergence  Rates  for  N  =  4,6,8 


Convergence  Rates 

N/Ne 

Ne  =  2  to  4 

Ne  =  4  to  8 

Ne  =  8  to  16 

N  =  4(co) 

3.2286 

4.1370 

4.1317 

3 

II 

£ 

3.2103 

6.3134 

6.8160 

N  =  6(co) 

5.8699 

6.0861 

6.0456 

3 

vO 

II 

£ 

7.3975 

8.9715 

7.9343 

N  =  8  (co) 

7.7170 

8.0165 

8.0291 

3 

00 

11 

£ 

8.9915 

10.1413 

9.8260 

Table  3.3:  Convergence  Rate  for  N  =  16 


Convergence  Rates 

N/Ne 

0 

CM 

II 

£ 

Ne  =  4  to  8 

N  =  16  (co) 

14.0437 

15.6555 

N  =  16  (p) 

16.5784 

17.5193 

When  testing  higher  order  polynomials,  we  can  increase  the  spatial  error  for  the 
given  function  and  decrease  the  time-step  in  order  to  ensure  that  we  can  neglect  any 
time-stepping  errors  during  testing  [1].  Doing  this  allowed  accurate  convergence 
rate  measurements  when  testing  16f/i  order  polynomials.  Figure  3.5  and  Table  3.3 
display  how  even  higher  order  polynomials  can  be  used  to  approximate  the  solution 
with  a  much  higher  convergence  rate.  However,  using  higher  order  polynomials 
comes  with  a  computational  cost  with  respect  to  time.  Higher  orders  are  much 
more  accurate,  as  seen  in  Figure  3.5,  but  the  computational  time  increases. 
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Figure  3.5:  Convergence  Rates  For  N  =  16  at  Ne  =  2,4,8 
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CHAPTER  4: 

Two-Dimensional  Discontinuous  Galerkin 


Moving  into  two  dimensions  is  obviously  more  challenging  than  one  dimension. 
The  first  portion  of  Chapter  4  focuses  on  an  example  two-dimensional  grid  with 
quadrilaterals  and  explains  the  formulation  of  the  metric  terms  needed  for  dis¬ 
cretization.  Section  4.2  is  the  discretization  of  (3.4),  (3.6),  and  (3.12)  in  two  dimen¬ 
sions.  As  discussed  in  Chapter  2,  working  with  integrals  and  mapping  the  physical 
domain  into  a  reference  domain  requires  a  change  of  variables.  In  two  dimensions, 
this  change  of  variables  produces  a  set  of  Jacobian  determinants  as  well  as  surface 
Jacobians  on  each  element.  Lastly,  we  apply  these  concepts  to  a  two-dimensional 
problem  on  a  washer  domain  with  curved  elements. 


4.1  Two-Dimensional  Grid 

Instead  of  a  one-dimensional  grid,  which  is  just  a  line  of  elements  in  Chapter  3, 
moving  into  two  dimensions  requires  a  two-dimensional  grid  of  elements.  Building 
the  grid  requires  a  change  of  variables,  from  the  x  and  y  spatial  coordinates  to  the 
£  and  r\  reference  coordinates  [9].  Figure  4.1  is  an  example  of  a  mapping  for 


Figure  4.1:  Example  Mapping  of  a  Quadrilateral  Element  in  Two  Dimensions  with  Degrees  of 
Freedom  for  N  =  2. 


a  quadrilateral  element  in  two  dimensions.  When  using  quadrilaterals  for  two 
dimensions,  we  have  ( N  + 1)2  LGL  points  in  both  the  Tj  and  £  directions,  as  depicted 
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by  the  red  X  symbols  in  Figure  4.1.  The  numbered  order  of  the  LGL  points  in 
Figure  4.1  is  important  for  the  storage  of  data.  For  example,  in  Section  4.2  the 
solution  is  stored  in  column  vectors 

<x)\ 

0)2 

co=  .  ,  p  = 

0)M 

where  there  are  M  =  (N  + 1)2  LGL  points  or  degrees  of  freedom. 

For  notational  purposes,  the  matrix  of  Jacobian  determinates  is  annotated  by  J  and 
the  matrix  of  surface  Jacobians  is  annotated  by  one  of  the  four  sides  of  the  element, 
such  as  Sp  for  side  one  of  the  element,  as  labeled  in  Figure  4.1.  Though  Figure  4.1 
is  an  example  of  elements  with  straight  boundaries,  we  develop  the  scheme  for 
general  curvilinear  quadrilateral  elements. 

4.1.1  Metric  Terms 

As  discussed  before,  we  need  to  conduct  a  change  of  variables  from  (x,  y)  to  (E,  ?/), 
where 


x  =  x(£,rj),  y  =  y{&,r\),  (4.1) 

and  we  assume  the  inverse  mapping  exists  with 

£  =  £(*,y),  *7  =  *?(*/ y)-  (4.2) 


Let  us  consider  a  single  quadrilateral  element,  such  as  the  one  in  Figure  4.1.  The 
following  explanation  for  the  metric  terms  is  based  on  Professor  F.X.  Giraldo's 
lecture  notes  [3].  Consider  the  differentials  defined  from  (4.1)  and  (4.2).  They  are 


dx  =  ^dE  +  ^-drr,  dE  =  ^-dx+^-dy, 
dE  or ]  dx  dy  “ 

dy  dy  dl]  dl] 

dy  =  -E-dE  +  -r-an,  dri=—dx+—dy, 
dE  dr j  dx  dy  “ 
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which  can  be  written  in  the  matrix  form 


dx 

dy 

dl 

dr\ 


"  dx 

dx 

dE 

di] 

dy 

dy 

dE 

di] 

~ dE 

dE' 

dx 

dy 

dr\ 

di] 

dx 

dy. 

dl 

dr\ 

dx 

dy 


(4.3) 

(4.4) 


The  Jacobian  and  the  inverse  Jacobian,  along  with  their  determinants,  from  (4.3) 
and  (4.4)  are 


J 


J 


-l 


dx  dx 
dE  di] 
dy  dy 

L  di  &}  j 
dE  dE 

dx  dy 
di J  di] 

L  dx  dy J 


.'-"'■SIM! 


.r*.«<rvS9-9S. 


dx  dy  dx  dy ' 


Focusing  on  (4.5),  taking  the  inverse  of  the  Jacobian,  J ,  yields 

(J)'1  = 


dy 


dx 
di] 

y  dx 


L  dE  dE 


_ i 

which  must  equal  JE  in  (4.6)  and  thus 


dE  dE 

dx  dy 
dr\  di] 
L  dx  dy 


dy 

di] 

dy 


dx 

di] 

dx 


L  dE  dE 


Equation  (4.8)  gives  the  metric  relations 

dy 


M  =  ]dd=_dV_  = 

dx  di]'  dx  dE,'  dy 


dx  di]  _  dx 
dr]'  Jdi/  =  dE' 


(4.5) 

(4.6) 


(4.7) 


(4.8) 


(4.9) 


which  arise  after  the  chain  rule  is  used  to  map  from  the  physical  domain  to  the 
reference  domain.  Lastly,  on  each  element,  we  store  the  Jacobian  determinants 
and  metric  terms  (4.9)  in  diagonal  matrices,  where  the  storage  order  is  as  the  LGL 
gridpoint  order  in  Figure  4.1. 
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4.2  Two-Dimensional  Discretization 

The  beauty  of  moving  into  two  dimensions  with  quadrilateral  elements  is  that  the 
constructions  of  the  mass  and  differentiation  matrices  are  relatively  simple.  The 
evaluation  of  the  surface  integral  over  each  element  requires  the  Kronecker  product 
of  the  one-dimensional  mass  matrix  with  itself.  This  allows  us  to  integrate  in  both 
the  £  and  rj  directions  for  each  element  and  is  annotated  by  M®M,  where  M  is  just 
the  one-dimensional  mass  matrix.  Additionally,  the  differentiation  matrices  are 

=  I0D, 

D;;  =  D<g>I, 

where  I  is  the  identity  matrix  of  size  (N  +  l)x(N  +  l)  and  D  is  the  one-dimensional 
differentiation  matrix. 

4.2.1  First  Two-Dimensional  Discretization  Equation 

The  discretization  of  (3.4)  is  very  similar  to  the  one-dimensional  discretization 
except  for  the  addition  of  the  Kronecker  product 

lT/(M®M)t-lT/(M®M)(y  =  0.  (4.10) 

clt 

Since  we  are  in  two  dimensions,  1  is  a  vector  of  ones  of  size  (N  + 1)2  and  is  needed 
because  in  (3.4)  we  are  integrating  against  the  function  ip(x)  =  1  [10]. 

4.2.2  Second  Two-Dimensional  Discretization  Equation 

The  discretization  for  (3.6)  requires  much  more  analysis  than  (3.4)  because  of  the 
gradients  and  surface  terms.  Due  to  this,  we  focus  on  the  left-hand  and  right-hand 
sides  separately.  Here  we  explicitly  introduce  the  differential  into  the  integral  for 
two  dimensions,  whereas  in  Section  3.4,  this  differential  was  implicit.  Focusing  on 
the  left-hand  side,  we  see  that 

XI  \  dr  r*+l 

+c2'Vip-'VpJdA  =  ipTJ(M<S>M)—  +  c2  J  J  [Vi/>  •  Vp] /d£drj. 
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where  Vi p  and  Vp  are  equal  to  the  following  after  the  change  of  variables: 


Vp 

dP*  dP  * 

dp 

dp*  dp  * 

dxl  +  dy^ 

dP 

dx  * +  dy] 

V:.  = 

dP*  dP  * 

dip 

dp*  dp  * 

dxt+  dy] 

dP 

li'  +  Ty] 

dp 
dp' 
dip 
dp  ’ 


Therefore,  Vip  •  Vp  yields  eight  separate  terms: 


Vi p-Vp 


dP  dp  dp  dp 
dx  dP  dx  dp  \ 


dP  dip  dp  dip 
dx  dP  dx  dp 


+ 


dP  dp  dp  dp 
dy  dd  dy  dp 


dP  dip  dp  dip 
dy  dP  dy  dp 


dP  dp  dP  dip  dP  dp  dp  dip  dp  dp  dP  dip  dP  dp  dp  dip 
dx  dP  dx  dP  dx  dP  dx  dp  dx  dp  dx  dP  dx  dp  dx  dp 


dP  dp  dP  dip  dP  dp  dp  dip  dp  dp  dP  dip  dP  dp  dp  dip 
dy  dP  dy  dP  dy  dP  dy  dp  dy  dp  dy  dp  dy  dp  dy  dp  ’ 


Focusing  on  the  first  term  of  the  second  equation  in  (4.11),  we  find 


r+1  r+1  (dp\dp (dp\dip  r+1  r+i/  dp\dp  t  dp\dip  1 

L  L  \i) L  <4-i2> 


■+1  r+1 


where  we  have  multiplied  in  JJ~X  to  place  a  Jacobian  determinate  with  both  of  the 
metric  terms.  By  substituting  in  the  metric  terms  from  (4.9)  and  discretizing,  we 
produce 


dy  dp  dy  dip 
dp  dP  dp  dP 


J~xdPdp 


'dy. 


dy. 


\dp 

=  V’X 


dp 


^rx(M0M)^ 

dp  dp  \ 


D&- 


(4.13) 


The  same  process  is  completed  for  each  of  the  eight  terms  in  (4.11)  to  find  the  final 
discrete  form  for  the  left-hand  side  of  (3.4)  to  be 


M)  ^  +c2ipTTp, 


(4.14) 
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where 


T={Dl 

^r1(M0M)^ 
dr]  dr] 

D,-D\ 

gr1(M®M)| 

~Dn 

||  /_1(M®M)| 

D,  +  Dj 

+DI 

jr-J1  (M0M) 
dr]  dr] 

D,-D\ 

~Dn 

—J  1  (M0M)  — 

De  +  D? 

(4.15) 


As  you  can  see,  the  discretization  for  two  dimensions  is  more  tedious  than  one 
dimension. 


The  right-hand  side  of  (3.6)  is  much  simpler  than  the  left-hand  side.  As  a  reminder, 
the  right-hand  side  of  (3.6)  is 

c2  (  rpn-  (' Vp)*ds . 

ddj 

Since  we  are  working  with  quadrilaterals,  we  integrate  the  four  boundaries  of 
each  element,  labeled  by  dD, j.  Figure  4.1  displays  how  the  four  boundaries  of 
each  element  are  numerically  labeled.  We  pull  the  information  from  each  side 
individually  by  using  the  following  operators: 

Li  =  Cg  <8)1,  L2  =  I® ejj,  L3  =  L4  =  I®eJ.  (4.16) 


The  four  L(  operators  (4.16)  isolate  the  correct  information  from  each  element 
boundary.  Additionally,  the  change  of  variables  for  the  right-hand  side  produces  a 
set  of  surface  Jacobians  found  in  either  the  £  or  ij  directions  by 


Sje  =  J 


Sj{  =  J 


£  =  1,3,  (4.17) 

£  =  2,4.  (4.18) 
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By  calculating  the  outward  unit  normal,  n,  for  the  sides  of  the  element  and  incor¬ 
porating  it  in  with  the  calculation  of  the  numerical  flux,  the  discretization  of  the 
right-hand  side  is 

<AT  [LjMSjj/p,  +  LT2MSj2fp2  +  LlMSpfp  +  LT4MSjifpi\,  (4.19) 

where  fp  =  n-  (Vp)*.  As  a  reminder  from  Section  4.1,  Sp,  Sn,  Sn,  and  Sy4  are  the 
matrix  of  surface  Jacobians  for  each  side  of  an  element,  as  labeled  in  Figure  4.1. 
This  same  notation  is  used  for  all  four  sides  of  the  flux  term,  /  ,  listed  in  (4.19). 
The  numerical  flux  for  two  dimensions  is  discussed  in  greater  detail  in  Section  4.3. 
Finally,  we  combine  the  two  sides  of  (3.6)  to  form  the  final  discretization  for  the 
second  variational  equation  to  be 

7(M®  M)  ^  +  c2Tp  =  c 2  (l[MS;1/pl  +  LT2MSj2fp2  +  LT3MSj3fp3  +  L[MS]Afp,) .  (4.20) 

As  a  reminder,  ip  does  not  appear  in  (4.20)  because  it  must  hold  for  all  ip,  as 
discussed  in  Section  3.3.1. 


4.2.3  Third  Two-Dimensional  Discretization  Equation 

The  discretization  process  of  (3.12)  is  similar  to  the  discretization  of  (3.6)  in  the 
previous  section.  Starting  with  the  left-hand  side. 


we  have  eight  terms  from  the  product  of  the  gradients 

dp  dip 
dx  dp 

d^  d  (dp  \d£dlp  dE  d  (dp  \  dp  dip 
dx  dp  \dt  )  dx  dE  dx  dp  \dt  J  dx  dp 

dE  d  (dp  \dEdlp  dE  d  (dp  \  dp  dip 

+  di/dE\dt~a,)di,dE+d^dE\dt~a,)di/^p 
dp  d  (dp  \  dE  dip  dE  d  (dp  \  dp  dip 

dy  dp  \dt  )  dy  dE  dy  dp  \dt  )  dy  dp 


Y7  I  vyldP  \_dE  d  (dp  \dEdlp  dE  d  dp 
^V1  dt  a>)  dxdEXdt  “IdxdE  dxdE\dt  a) 


(4.21) 
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Focusing  on  the  first  term  of  (4.21),  we  find  the  discrete  form  to  be 


n+1 
1 


xpTD\ 


d4I-4  <4-22> 


Expanding  the  discretization  for  all  eight  terms  in  (4.21)  yields  the  full  discrete 
form  for  the  left-hand  side  of  (3.12), 


T  ( dp  T  dp  T 

xbTT  H--co  =  \bT T -j- -  \bT T co , 
r  \dt  r  dt  r 


(4.23) 


where  T  is  equal  to  Equation  (4.15). 

Let  us  focus  on  the  right-hand  side  of  (3.12), 


(  (Vxp-n)(co*  -co)ds, 

J  <9Q, 


(4.24) 


and  discretize  side  one  of  a  single  element  (see  Figure  4.1).  Conducting  a  change 
of  variables  yields 


J'  Vip-n  fC0Sjdc;r 

where  the  fa,  =  ( co *  -  co).  For  notational  purposes,  allow 

Dx 


(4.25) 


D 


y 


dl  ^  dl] 
d^D^+d^D'\ 
dl  dll 

[  TyD^TyD\ 


The  discretization  of  (4.24)  for  side  one  becomes 

\T 


r 


(L\Dx)  nx  +  ^LiDy^j  riy  MSjif ay 
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where  the  flux,  f  y  is  from  side  one  of  an  element  as  annotated  in  Figure  4.1.  A 
similar  calculation  for  the  remaining  sides  gives 


<PT  (dJ  LjVr - + + 4>r  (dJl  + D  Jb,) 

+4’r  [DlLln  „  +  Djljny)i MS, 3/^  +  t/>T  (dJl[hv  +  D  MS^. 

Equation  (4.26)  is  the  discrete  form  for  the  right  side  of  (3.12).  Combining  (4.26) 
and  (4.23)  produces  the  final  discrete  approximation  for  the  third  variational  equa¬ 
tion  to  be 


T(o  —  (DxL1nx  +  DyL 1  ny^MS j\f  a)1  + 

+  (DjL \nx  +  DTyLlny)MS^f(o3  +  ( 


DTxLlnx  +  DTyLlny 
DTxLlnx  +  DTyLlny 


)MSq/w4. 


(4.27) 


4.2.4  Combination  of  Two-Dimensional  Discretization  Equations 

As  in  one  dimension,  we  have  three  discrete  equations  and  our  goal  is  to  get  the 
problem  into  the  form  of  two  equations  and  two  vector  unknowns.  By  combin¬ 
ing  (4.10)  and  (4.27),  we  find 

(r+l/(M®M))^  =  (r+  )a>  +  +  Djt[«y)MS)1/0,1 

+  (DjLf«,  +  Djq'H,)MS/2/aa  +  (Djl|»I  +  DjlJ»y)MS;3/<03  <4'28> 

+  (dJlJm*  +  DjJLjny)MSj4fa)y 

where  1  is  a  matrix  of  ones  since  (4.10)  produces  a  scalar.  Additionally,  for  the  sec¬ 
ond  equation,  (4.20)  can  be  re-arranged  to  isolate  %  for  the  final  discrete  equation 
of 

7(M®M)  A  =  ((LfMS(1/fl  +  LT2MSj2fp2  +  L2MSjsfp3  +  LptfS, 4/p4)  -  Tp)  .(4.29) 

Equations  (4.28)  and  (4.29)  are  the  two-dimensional  ordinary  differential  equations 
that  we  solve  numerically  using  the  RK54  scheme  [11]. 
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4.3  Two-Dimensional  Numerical  Flux 

The  numerical  flux  for  two  dimensions  is  conceptually  similar  to  the  numerical  flux 
for  one  dimension,  except  in  two  dimensions  we  have  to  account  for  more  degrees 
of  freedom,  boundaries,  and  the  element  normals.  In  Section  3.5,  we  accounted 
for  the  discontinuity  that  existed  between  the  two  different  approximations  at  the 
boundary  point.  In  this  section,  we  account  for  the  discontinuity  that  exists  between 
two  different  approximations  at  multiple  boundary  points  along  the  surface  of  two 
neighboring  elements.  Thus,  the  numerical  flux  in  two  dimensions  is  an  extension 
of  the  numerical  flux  from  Section  3.5.  Figure  4.2  is  an  example  depiction  of  a 
two-dimensional  flux  boundary  for  a  quadrilateral  element. 
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Figure  4.2:  Example  Two-Dimensional  Flux  Boundary  on  a  Quadrilateral  with  N  =  2. 


In  two  dimensions,  the  element  to  be  considered  is  the  (-)  side  and  the  neighboring 
element  is  the  (+)  side.  This  is  the  same  for  the  (-)  representing  the  element  of 
focus's  gridpoints  and  the  (+)  representing  the  neighboring  elements  gridpoints. 
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In  referencing  Figure  4.2,  the  center  element  requires  a  flux  to  be  computed  with 
its  neighbors  at  the  boundary  LGL  points.  Once  again.  Figure  4.2  is  an  example 
for  discussion  purposes;  the  actual  elements  are  curvilinear  quadrilaterals  and  the 
element  normals  at  the  LGL  points  are  required  in  the  numerical  flux  computation 
of  (ft  ■  Vp)*.  The  formulation  of  the  form  of  the  central  flux  and  upwind  flux  can  be 
found  in  Section  3.5.  The  central  flux  equations  for  two  dimensions  are 

co*  =  ^  (co~  +co+),  (4.30) 

n  ■  (Vp)’  =  \  (n~  ■  (Vp~)  -  n*  ■  (Vp*)) ,  (4.31) 

where  (4.30)  element  normals  are  represented  in  (4.27).  The  element  normals  in  two 
dimensions  are  outward  of  each  element.  For  example,  n~  is  the  normal  directed 
away  from  the  (-)  boundary  toward  the  neighboring  (+)  element  and  n+  is  directed 
toward  the  element  of  focus  (i.e.,  the  center  element  in  Figure  4.2).  Typically,  the 
central  flux  is  an  average  of  both  sides  of  an  element  boundary,  but  the  reader  may 
notice  that  (4.31)  has  a  minus  instead  of  a  plus  sign.  The  reason  for  this  is  that 
ft+  =  -ft-,  and  thus  the  minus  sign  actually  produces  an  average.  As  a  reminder, 
the  numerical  flux  is  incorporated  in  the  discretization  equations  by  fp  =  n-  (Vp)* 
from  Section  4.2.2  and  fw  =  (co*  -  co)  from  Section  4.2.3.  Similar  to  one  dimension, 
the  central  flux  is  easy  to  understand  and  produces  a  stable  algorithm,  but  it  does 
not  take  into  account  the  physical  propagation  of  information  like  the  upwind  flux. 

The  upwind  flux  equations  for  two  dimensions  are 

<*<’  =  ( (to~  +  a>*)  -  J  (n*  ■  (Vp*)  +  n~  ■  (Vp-)) ,  (4.32) 

« ' (Vp)'  =!(«-■  (Vp“) - «+  -  (Vp+))  +  1  ( co* -a,').  (4.33) 

Equations  (4.32)  and  (4.33)  are  formulated  by  the  decoupling  of  the  wave  equation 
into  two  one-way  advection  equations  at  the  boundary  points  through  the  charac¬ 
teristic  variables.  This  is  the  same  concept  from  Section  3.5,  but  in  two  dimensions, 
there  are  three  characteristic  variables.  However,  the  third  characteristic  variable 
is  associated  with  a  zero  eigenvalue  and  does  not  propagate  information  across 
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element  boundaries.  As  depicted  in  Figure  4.2,  r1  and  are  example  characteristic 
variables  propagating  information  across  one  elemental  boundary. 


4.4  Two-Dimensional  Discontinuous  Galerkin  Results 

As  a  test  problem  we  consider  a  washer  domain  with  curved  elements  as  shown 
in  Figure  4.3.  Along  the  inner  and  outer  radius  of  the  washer,  we  impose  the  free 
surface  boundary  conditions  of  p  =  0.  We  use  the  exact  solution  of 


dp 

Tt 


p(x, y,  t)  =  sin  ( nt  -  / 66)  Jp  (r(x,  y)) , 
co(x,  y,  t)  =  n  cos  (nt-  fid)  Jp  (r(x,  y)) , 


with  the  conversion  equations 


Kx,y)  =  J- 


x2  +  y2, 


6(x,y )  =  cos 


-i| 


x 


j(x,y), 


(4.34) 

(4.35) 


for  mapping  between  polar  and  Cartesian  coordinates.  Here,  Jp  is  a  Bessel  function 
of  the  first  kind  with  parameter  /3  and  n  is  an  integer;  we  use  /3  =  4  and  n  =  1.  To 
ensure  (4.34)  is  equal  to  zero  at  the  inner  and  out  radius  for  all  time,  we  enforced  the 
inner  and  outer  radius  of  the  washer  grid  to  be  the  second  and  fourth  roots  of  the 
Bessel  function.  For  interested  readers,  the  DG  implementation  for  two  dimensions 
can  be  found  in  Appendix  C. 


We  initially  tested  the  implementation  on  a  square  grid  of  quadrilateral  elements 
with  periodic  boundary  conditions.  After  the  successful  implementation  on  a 
square  grid,  we  built  and  tested  the  algorithm  on  a  washer  grid  with  curved 
elements  using  the  exact  solution  of  (4.34)  and  (4.35)  evaluated  at  t  =  0  as  the  initial 
condition.  To  enforce  the  free  surface  boundary  condition  of  p  =  0,  we  used 


C 0+  =  -cv  ,  (n+-(Vp+))  =  -(n  -(Vp  )) 


in  the  numerical  flux  routine  for  the  points  along  the  inner  and  outer  radius  of  the 


44 


2 


washer.  Once  again,  a  visual  depiction  of  an  example  washer  grid  with  Ne  =  8 
semi-curved  elements  is  seen  in  Figure  4.3.  Initially,  testing  the  entire  washer  grid 
became  too  computationally  expensive  with  the  large  number  of  elements.  To 
counter  this  problem  and  still  maintain  integrity  of  the  results,  we  evaluated  one 
quarter  of  the  washer  by  enforcing  j  periodic  boundary  conditions  for  the  left  and 
right  boundaries  through  the  bessel  function.  By  letting  /3  be  any  multiple  of  four, 
we  enforced  the  j  periodic  conditions.  For  example,  in  referencing  Figure  4.3,  this 
means  boundaries  1  and  3  from  elements  1  and  2  are  equal. 

Table  4.1  lists  the  different  polynomial  orders  and  elements  tested  using  the  initial 
condition  with  (4.28)  and  (4.29)  and  evaluated  through  a  RK54  iterative  method. 
Figure  4.4  shows  the  log  of  the  error  vs.  log  of  the  number  of  elements  for  N  = 
2,  4,  6,  8  calculated  using  the  global  L2  error  norms.  As  expected  for  both  co 
and  p,  Figure  4.4  displays  increasing  convergence  rates  as  the  polynomial  order 
increases.  Once  again,  the  higher  the  order  of  the  local  approximation,  the  faster 
the  convergence  rates  are  due  to  the  the  global  error  being  dependent  on  the 
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Table  4.1:  Discontinuous  Galerkin  Tested  Information 


Polynomial  Orders  (N) 

Number  of  Elements  (Ne) 

N  =  2, 4,  6,  8 

Ne  =  4, 16,  64, 256 

N  =  10, 12 

Ne  =  4, 16,  64 

Convergence  Plot  for  w 


Figure  4.4:  Convergence  Rates  For  N  =  2,  4,  6,  8 

polynomial  order  [1]  (as  discussed  in  Chapter  3  with  (3.32)). 

Table  4.2  displays  the  convergence  rates  corresponding  to  Figure  4.4.  As  seen, 
the  convergence  rates  increase  with  increasing  polynomial  order.  Using  the  same 
space  of  functions  for  co  and  p,  the  two-dimensional  results  yielded  convergence 
rates  being  near  its  associated  polynomial  order,  with  co  being  near  N  and  p  being 
N  +  \  or  N  + 1.  Once  again,  most  DG  methods  are  expected  to  be  of  the  order  N, 
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Table  4.2:  Convergence  Rates  for  N  =  2, 4,  6,  8 


Convergence  Rates 

N/Ne 

Ne  =  4  to  16 

Ne  =  16  to  64 

Ne  =  64  to  256 

N  =  2(o)) 

1.9240 

2.0570 

2.1274 

.a, 

CM 

II 

£ 

1.9663 

2.9243 

2.8680 

N  =  4(oj) 

3.3438 

4.2448 

4.0822 

.a. 

II 

£ 

4.2739 

5.1808 

4.4615 

N  =  6  (o)) 

5.2521 

6.0193 

6.0944 

.a. 

V© 

II 

£ 

6.6670 

6.7889 

6.4519 

N  =  8(co) 

6.9754 

7.9190 

8.0543 

.a. 

00 

II 

£ 

8.5157 

8.7873 

8.5069 

N+  orN  +  1  [1].  These  results  are  different  from  the  one-dimensional  convergence 
rates,  where  co  is  N  and  p  is  N  +  2.  Understanding  this  difference  requires  further 
analysis  and  is  an  area  for  future  research. 


Table  4.3:  Convergence  Rate  for  N  =  10, 12 


Convergence  Rates 

N/Ne 

Ne  =  4  to  16 

Ne  =  16  to  64 

N  =  10  (co) 

8.7628 

10.0287 

N  -10  (p) 

10.2080 

10.6994 

N=  12  (co) 

10.6378 

12.0772 

N  -12  (p) 

12.0461 

11.8315 

To  use  higher-order  polynomials,  we  decreased  the  number  of  elements  in  order 
to  avoid  reaching  machine  precision  with  the  first  datapoint.  Figure  4.5  shows 
the  error  and  Table  4.3  gives  the  convergence  rates  for  N  —  10  and  N  =  12.  The 
convergence  rates  are  still  on  the  order  of  the  polynomial,  but  the  N  =  12  case 
seems  to  be  reaching  machine  precision.  Both  Figure  4.5  and  Table  4.3  display  how 
higher  order  polynomials  can  be  used  to  approximate  the  solution  with  a  much 
higher  convergence  rate. 

Using  higher  order  polynomials,  especially  on  complex  geometries  with  curved 
elements,  comes  with  a  computational  cost  with  respect  to  time  for  approximating 
the  solution.  Higher  orders  are  more  accurate,  as  seen  in  Figure  4.5,  but  the 
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Figure  4.5:  Convergence  Rates  For  N  =  10, 12 

computational  time  increases  since  the  matrices  to  be  inverted  are  larger  with 
higher  polynomial  orders  and  require  a  smaller  time-step.  In  the  code  we  use  a 
time-step  of  At  ~  A. 
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CHAPTER  5: 
Energy  Conservation 


Proving  stability  of  a  method  can  be  accomplished  through  an  energy  analysis. 
The  following  energy  analysis  is  similar  to  the  method  employed  by  Appelo  and 
Hagstrom  [2].  By  conducting  the  energy  analysis  on  an  element  with  neighboring 
elements,  we  are  able  to  find  out  if  energy  is  dissipated  or  conserved  throughout 
the  system.  If  energy  is  dissipated,  then  we  know  that  all  of  the  eigenvalues  of 
the  system  have  a  negative  real  part  and  thus  the  ODE  is  stable.  If  energy  (E)  is 
conserved,  then  ^  =  0  and  all  the  eigenvalues  are  purely  imaginary.  If  E k  is  the 
energy  on  an  element  k,  then  the  total  energy  is 

Ne 

E  =  Y4Ek,  (5.1) 

k= 1 

and  we  want  to  show 


dE  dEk 
— = >  — ^<0 
dt  '  ^  dt 

k=  1 


(5.2) 


in  order  to  prove  that  energy  is  dissipated  throughout  the  system  as  time  progresses. 


5.1  Basic  Theory  of  Energy  Conservation 

Let  us  look  at  the  continuous  energy  equation  on  the  domain  x  G  [a,  b]  for  one 
dimension: 


(5.3) 


with  pt  =  co,  cot  =  c2pxx,  and  c2  =  We  assume  c.  A,  and  p  are  constants  on  the 

domain.  As  a  reminder,  for  clarification  purposes,  %  —  pt  for  a  short  notation  form 
and  this  follows  through  for  the  remaining  terms  in  Section  5.1.  The  first  portion 
of  (5.3)  is  the  kinetic  energy  and  the  second  portion  is  the  potential  energy  for  the 
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system.  By  taking  the  time  derivative  of  (5.3) 


^  f'l  1 

Et  =  I  jCOtCO  +  -pxtPxdx, 


(5.4) 


we  can  now  begin  to  manipulate  (5.4)  to  show  (5.2).  Specifically,  we  want  to 
manipulate  (5.4)  in  order  to  select  our  boundary  conditions  to  ensure  Ef  <  0.  Let  us 
substitute  ay  =  c2pxx  and  co  =  pt  into  (5.4) 


Ef  = 


1  2  1 
jcyxxpt  +  -pxtpxdx. 


(5.5) 


With  c2  =  jj,  it  follows  that  (5.5)  becomes 

1  rb 

Et  -  —  Pxxpt+ Pxtpxdx,  (5.6) 

P  Ja 

and  by  conducting  integration  by  parts  on  the  last  portion  of  (5.6)  we  find  the 
following: 


Et-~ 

P 


Jrb  r*b 

PxxPtdx  +  pxPt  \„-  I  VxxVtdx 
a  Ja 


(5.7) 


The  two  remaining  integrals  cancel  and  we  have  that 

1  1 

Ef  =  ~pxpt  l!,=  -  [pXCO  |f,  -pxCO  | a], 


(5.8) 


where  we  can  now  choose  the  boundary  conditions  to  ensure  Ef  <  0.  For  example, 
with  periodic  boundary  conditions,  it  should  be  obvious  that  (5.8)  equals  zero  when 
a  —  b  and  energy  is  conserved. 


5.2  Two-Dimensional  Energy  Analysis 

Section  5.1  is  a  relatively  simple  example  of  the  analysis  of  energy  for  the  continu¬ 
ous  one-dimensional  system.  Moving  into  two  dimensions  requires  more  tedious 
analysis  and  also  depends  on  the  conditions  set  within  the  system.  For  example. 
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when  using  an  upwind  flux,  we  want  to  show  that  energy  is  dissipated  and  when 
using  a  central  flux,  we  want  to  show  that  energy  is  conserved.  Let  us  initially 
investigate  a  general  two-dimensional  case  before  moving  into  the  flux  analysis. 
The  energy  equation  for  two  dimensions  is 

Ek  =  LM  +  Tp{pi+^dA  (5'9> 

where 


dp  dE  dp  dl]  dp 
dx  dx  dE  dx  dp ' 


and 


/  dp  \2  dp  /  dE  \2  dp  dp  dE  dp  dp  dp  dp  dE  dp  dp  /  dp  \2  dp 
[dxj  dE\dx)  dE  dEdxdxdp  dp  dx  dx  dE  dp\dx)  dp' 


Moreover, 


is  formulated  through  the  same  method  as 


respect  to  y.  The  change  of  variables  of  (5.9)  yields 


*-rrw 


•% 

dp 

%! 


dE  dp  dE  dp 


dp 

dE 1 

L 

dp  dp 


<d_e 

i  dx , 


21 


+ 


dx  dx  dy  dy  \  dp  dp 
dEdp 


<d_E 
\dlJ 
dE  dp 


dp 

dE 

dp 


dAdJl  + 

dx  dx  dydy\dE 


(dZ)2  ,| 

'dE\2 

dp 

\dx)  1 

.dy)  _ 

dp 

(5.10) 


and  the  discretization  of  (5.10)  is 


2  A 


2  p' 


(5.11) 
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where  T  is  the  same  simplifying  variable  used  in  Equation  (4.14)  from  Section  4.2.2. 
Taking  the  time  derivative  of  (5.11)  we  have 


dEfc 

dt 


1  Ida)'' 
2A\  dt  i 


+ 


2 p  \  dt 


J  (M®M)co+  —coJJ(M®M) 
2  A 


Fp  +  2pPTT\ 


which  can  be  simplified  to 


dEk 

dt 


Lt/(M®M) 

A 


T 


(5.12) 


Upon  inspection  of  (5.12),  we  can  substitute  (4.27)  and  (4.29)  in  for  J(M®M)|jj 
and  from  Chapter  4  to  produce 

jp  2 

-df  =  Ja)T  liLlMSnfpi  +  L\ lMSj2fp2  +  LT3MSj3fp  3  +  L4TMS/4/p4)  -  Tp] 

+  -ppJ[Tcv  +  (DjL[«x  +  Djl[«y)MS;1/a,1  +  (DTxLT2nx  +  DTyLT2ny)MSj2fJ5-13) 
+  (DIL3«t  +  DyL^ny^MS  j3f  co3  +  (dJlJm*  +  D[Ljn1/)MS/4/fi)4] . 

For  notational  purposes,  let 

Fp  =  (%MSnfpl  +  LT2MSj2fp2  +  LT3MSj3fp3  +  L4TMS;4/p4) , 

F(l,  =  (DlLlnx  +  D^Llny)MSj1fwl  +  {pl.Llnx  +  Dj/Llny)MSj2fw  2 

+  ( D^L3nx  +  DyL3ny^MS  j3f w3  +  ( D^L^nx  +  D^L^riyj  MSj^f aj4. 

2  -i 

Therefore,  Equation  (5.13),  with  the  substitution  of  j-  -  ±  becomes 

=  ^c°T(fp~tp)  +  ^PT^Tco  +  f^'  (5-l4) 
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and  since  coTTp  =  pTFco,  (5.14)  is  further  simplified  to 


f  =  I(o .WrB).  («5) 

As  you  can  see,  (5.15)  is  dependent  on  the  flux  and  boundary  information  from 
each  element,  as  included  in  Fp  and  Fa).  Knowing  that 

Cl)  g  -  Cl)  L-l  £  f 

(ne  •  (Vp^))T  =  pT  (DlLT£nx  +  Djljny) , 


for  €  =  1,2,3, 4,  (5.15)  becomes 


dE  1 

-Jf  =  -  [w[MS/l/pl  +  a)2MSl2fp2  +  C°3MSpfp3  +  a)lMSj^fpA 

+  («1  •  (Vp1))TMS/1/a,1  +  («2  •  (Vp2)f  MS  j2fa)2 
+  («3  •  (Vp3))TMS/3/&, 3  +  («4  •  (Vp4))TMS74/al4  • 


Recall  that  fp-n  ■  (Vp)*  and  fa)  =  ( co *  -  co)  from  Section  4.2.2  and  Section  4.2.3,  thus 
we  can  consolidate  the  boundary  information  for  each  side  of  an  element 


dEp 

dt 


(o^MSp  (n\  •  (Vp4))  +(«r(Vp1))  MSp  (co*  -  co\) 

+  ( 'coT2MSj2  («2  •  (Vp2))  +  (n2  ■  (yp2))T MS j2  (co*  -  co2) 
+  (a)lMSj3(n3-(yp3))  +  («3  •  (Vp3))T MS j3(co*-co3) 
+  ( o)4MS;-4  (n4  •  (Vp4))  +  («4  •  (Vp4))TMSj4  (o>*  -  co4) 


(5.17) 


Equation  (5.17)  is  in  the  final  form  that  we  use  to  conduct  the  stability  analysis. 
Summing  (5.17)  for  all  elements,  we  can  analyze  the  energy  for  the  whole  domain. 
In  what  follows,  we  prove  stability  with  the  boundary  conditions  imposed  from 
Chapter  4  for  the  central  flux  and  the  upwind  flux. 
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5.2.1  Central  Flux  Analysis 

In  order  to  simplify  the  analysis,  we  focus  on  a  single  side  of  an  element  with  its 
adjacent  neighboring  element.  This  concept  is  depicted  in  Figure  3.2  in  Chapter  3. 
For  notation  purposes,  the  (-)  side  is  the  element  of  focus  and  the  (+)  side  is  the 
neighboring  element.  Additionally,  co~  represents  the  solution  at  the  degrees  of 
freedom  on  the  associated  side  of  the  quadrilateral  element  as  depicted  in  Figure  4.2. 
This  same  notation  applies  for  (5.19)  and  (5.20).  The  equations  for  the  central  flux 
for  a  generic  side  without  boundary  conditions  are 


+10+), 

(5.18) 

(Vp)=l[«--(Vp-)-n+-(Vp+)], 

(5.19) 

(Vp,)  =  l[«+-(Vpt)-«--(Vp-)]. 

(5.20) 

Since  we  are  working  with  one  side  of  a  quadrilateral  element,  we  add  together 
the  terms  from  one  edge  of  the  element  with  the  terms  from  the  neighbor's  corre¬ 
sponding  edge 


I si- l—i -I  p  rp  rp  1 

=  [(^+)  MS;-  ( n *  ■  (Vp)*)  +  (n+  ■  (Vp+))TMS;  («•  -  w+)]  (g  2J) 

+  [{co~)T MS j  (n~  •  (Vp)*)  +  («-  •  ( Vp~))TMSj  (m*  -  oT)\ . 

Substituting  in  (5.18),  (5.19),  and  (5.20)  for  the  flux  portions  in  (5.21)  yields 
dE±  i 

[(co+f  MS  j(n+  •  (Vp+))  -  (n+  •  (Vp+))T  MS jCO+  +  (n+  ■  (Vp+))T  MS jCO~ 

-  (m“)T MS j (n+  ■  (Vp+))  +  (m“)TMS;- {n  ■  (Vp“))  -  {n  ■  (Vp“))T MS jCO~  (5’22) 
+  {n  ■  (Vp“))T  MS j(o+  -  ((o+)T  MS  j  («"  •  (Vp“))] . 

dEf 

Inspecting  (5.22),  all  of  the  terms  cancel  out  to  show  -A  =  0  for  one  element 
boundary.  This  can  be  expanded  to  all  four  boundaries  of  the  quadrilateral.  Thus, 
energy  is  conserved  and  the  ODE  is  stable  when  using  the  central  flux. 
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5.2.2  Upwind  Flux  Analysis 

As  discussed  at  the  end  of  Section  3.5,  the  upwind  flux  consists  of  the  central  flux 

plus  an  upwinding  portion.  This  is  evident  in  (3.28)  and  (3.29).  Therefore,  since  we 

dEf 

know  that  energy  is  conserved  (i.e.,  -A  =  0)  for  the  central  flux  portion,  we  only 
investigate  the  upwinding  portion.  The  equations  for  the  upwinding  section  are 


K  =  - \  [n+  ■  (Vp+)  +  n  ■  (Vp-)] , 

(5.23) 

(»+,(Vp)*)M  =  ^(oT-<y+), 

(5.24) 

(«"-(Vp)*)„  =  ^(oA-oT), 

(5.25) 

which  is  substituted  into  the  elemental  boundary  equation  (5.26)  for  the  upwinding 
portion. 


nr)  =  [(‘*’+>Tms<  •  W)» + (”+  '  (Vp+))TMS,a>;] 

+ [(«-)tms,  („-  ■  (Vp)*)„ + («-  ■  (Vp-))TMS,ffl;,] . 

The  substitution  of  (5.23),  (5.24),  and  (5.25)  into  (5.26)  yields 


dE? 

_ k_ 

dt 


T 

=  —  [(<u+)T MS j(co~  -  co+)  +  (co~)T MS j(co+  -  <u-)] 

- 1  [(«+  •  (Vp+))T MS,  [ n *  ■  (Vp+)  +  n  ■  (Vp-)] 

+  («-  •  (Vp-))T MS j  [n+  •  (Vp+)  +  n  •  (V/T)]] . 


(5.26) 


(5.27) 


Factoring  out  a  negative  from  the  first  portion  of  (5.27)  makes  (to  -  co+)  become 
( co+  -  co~).  Therefore,  (5.27)  can  be  re-written  to 


Yc[(co+-co  )T  MS  j(to+  -  to  )] 

\  [(n+  •  (Vp+)  +  n  •  (Vp~))T MS  j  (n+  •  (Vp+)  +  n  •  (Vp"))] . 


(5.28) 
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Inspection  of  (5.28)  displays  two  terms: 


[(<d+-<d  )T MS j{co+  - co  )], 

[(n+  •  (Vp+)  +  n  ■  (Vp~))T MS j  ( n+  •  (Vp+)  +  n  •  (V/T))] 

that  are  both  positive  definite.  MSj  in  both  sections  is  always  positive;  therefore,  re¬ 
gardless  if  ( co+  -  co~)  or  (n+  •  (Vp+)  +  n~  •  (Vp-))  produce  negative  or  positive  values, 
the  square  of  each  is  always  be  positive.  Furthermore,  with  the  negative  in  front 
of  both  sections  of  (5.28),  <  0.  The  combination  of  with  ^  from  the 

central  flux  section  (Section  5.2.1)  yields  a  negative  value  per  face.  The  summation 
of  all  the  elements  in  the  domain  shows  Equation  (5.2)  to  be  true.  Therefore,  this 
method  is  stable  for  both  the  central  flux  and  upwind  flux  on  complex  geometries 
with  curved  elements. 

5.2.3  Boundary  Condition  Analysis 

In  Chapter  4,  we  tested  the  method  on  a  washer  mesh  with  curved  elements  and 
the  boundary  conditions  of  p  =  0  on  the  inner  and  outer  radius  of  the  washer.  In 
order  to  implement  the  boundary  conditions,  the  following  constraints  were  set  at 
the  inner  and  outer  boundaries  of  the  washer 

co+  =  -co~r  (n+-(Vp+))  =  -(fT-(Vp-)).  (5.29) 

Additionally,  we  are  only  concerned  with  the  (-)  portion  of  (5.21)  for  the  central 
flux  and  the  (-)  portion  of  (5.26)  for  the  upwind  flux  since  there  are  no  neighboring 
elements  ((+)  portions)  at  the  inner  and  outer  radius  of  the  washer. 


Boundary  Conditions  with  Central  Flux 

Once  again,  substituting  in  the  central  flux  Equations  (5.18)  and  (5.19)  for  the  (-) 
portion  of  (5.21)  yields 

=  (a r)TMS;(l  [«-  ■  (Vp-)  -  «+  ■  (Vp+)])  +  («-  ■  { Vp-))TMS;( \ (w*  -  or)). 
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which  becomes 


[2(„-.(Vp-))])  +  („-.(Vp-))TMS, 

after  the  substitution  of  the  boundary  conditions  listed  in  (5.29).  Simplifying  (5.30) 
produces 

dlij~ 

=  (co~)T MS j(n~  ■  (Vp-))  -  {n~  ■  (Wp~))T MS j(co~)  =  0,  (5.31) 

and  energy  is  conserved  with  the  boundary  conditions  of  p  =  0  using  a  central  flux. 

Boundary  Conditions  with  Upwind  Flux 

For  the  upwind  flux  boundary  condition  analysis,  we  use  the  same  process  as 
above,  except  we  use  (5.23)  and  (5.25)  for  substitution  into  the  (-)  portion  of  (5.26). 
This  yields 


dE, 

_ k 

dt 


=  (co-)TMSi[  i 


( dE1 

_ k_ 

i  dt 


(a)  )TMSj^(co+-co  ))  +  («  (Vp  ))TMS;(-^[n+-(Vp+)  +  n  -(Vp  )]) 


which  becomes 


( dE1 

_ k_ 

i  dt 


(co-f  MS,(±-c (-20.-))  +  («-  •  (Vp-))TMS7-(-  C- [-«-  •  (Vp-)  +  n  ■  (Vp“)]) 


and  is  simplified  to 


j  =  -- [(a)-)TMSjO)-]  <  0.  (5.32) 

/  u  c 

Once  again,  the  combination  of  (5.31)  and  (5.32)  produces  a  negative  result  for  the 
upwind  flux  boundary  condition  analysis. 

The  global  energy  dissipation  rate  is  the  sum  of  the  energy  dissipation  rates  at  all  of 
the  faces.  Therefore,  for  both  the  central  and  upwind  fluxes,  the  method  discussed 
in  Chapter  4  is  stable. 
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CHAPTER  6: 
Conclusion 


The  complexities  of  solving  PDEs  through  Finite  Difference  Methods,  Finite  Volume 
Methods,  and  Finite  Element  Methods  is  and  will  continue  to  be  an  area  of  ongoing 
research.  In  this  paper,  we  explored  a  new  Discontinuous  Galerkin  method  for 
approximating  a  second  order  wave  equation  with  curved  elements  in  complex 
geometries.  The  beginning  two  chapters  discuss  the  tools  needed  for  constructing 
the  DG  method,  while  Chapter  3  and  Chapter  4  tested  the  method  in  one  and  two 
dimensions.  Chapter  5  proved  energy  stability  of  the  method  through  an  energy 
analysis. 

In  Chapter  3,  convergence  rates  on  a  one-dimensional  grid  for  both  low  and  high 
polynomial  orders  yielded  rates  near  the  associated  polynomial  order.  The  results 
yielded  convergence  rates  with  co  being  on  the  order  of  N  and  p  being  on  the  order 
of  N  +  2,  which  is  higher  than  most  DG  methods  that  yield  results  on  the  order  of 
N,  N  +  and  N  +  1.  Additionally,  using  higher  polynomial  orders  required  fewer 
elements  for  a  given  error  level. 

In  two  dimensions,  we  tested  the  problem  on  a  washer  grid  with  curved  elements. 
In  this  case,  the  inner  and  outer  radius  had  the  free  surface  boundary  condition 
of  p  =  0  implemented  through  the  numerical  flux  routine.  The  two-dimensional 
results  yielded  convergence  rates  different  than  the  one-dimensional  method  with 
co  and  p  being  near  its  associated  polynomial  order  N  or  N  +  Similar  to  the  one¬ 
dimensional  results,  testing  high  polynomial  orders  achieved  machine  precision  at 
a  faster  rate  with  less  elements. 

Chapter  5  explored  the  energy  dissipation  of  the  method,  with  and  without  the 
imposed  boundary  conditions  of  p  =  0,  for  the  central  flux  and  upwind  flux.  In  all 
cases,  the  method  either  conserved  energy  (^  =  0)  or  dissipated  energy  (  jk  <  o). 
Therefore,  the  eigenvalues  of  the  system  have  a  negative  real  part  proving  that  the 
ODE  is  stable  and  the  method  is  viable. 
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6.1  Future  Work 

Future  research  into  this  Discontinuous  Galerkin  method  can  be  conducted  on  a 
variety  of  topics.  These  topics  include,  but  are  not  limited  to,  the  following: 

•  Accuracy:  As  discussed  in  Chapter  3  and  Chapter  4,  the  convergence  rates 
listed  are  observations;  proving  these  results  could  lead  to  further  insight  and 
improvement. 

•  Efficiency:  Is  it  possible  to  template  the  curvilinear  elements  so  that  the  mass 
matrix  is  the  same  for  all  elements  by  modifying  the  approximation  space  for 
each  element?  The  benefit  of  doing  this  is  the  ability  to  store  one  mass  matrix 
for  all  elements  [12, 13]. 

•  Coupling:  Development  of  the  numerical  coupling  procedures  for  elastic- 
acoustic  interfaces  within  this  second  order  form  [14]. 

•  Form:  Relationship  between  this  second  order  formulation  and  the  standard 
first  order  formulation  of  the  acoustic  wave  equation. 

•  Extension:  Consider  this  method  for  other  equations,  such  as  the  Einstein 
equations  governing  black  hole  dynamics  where  there  are  10  equations  in  the 
second  order  form  versus  50  equations  in  the  first  order  form  [15]. 
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APPENDIX  A: 
Interpolation  and  Integration 


The  following  are  the  main  MATLAB  codes  for  Interpolation  and  Integration  from 
Chapter  2. 


A.l  Interpolation 


1 


3 


5 


7 


9 


11 


13 


15 


17 


19 


21 


23 


25 


27 


o/  o/ 

/o - /o 

%This  is  the  main  driver  for  Interpolation  . 

%Written  by  Benjamin  Davis 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%Synopsis :  Conducting  Interpolation  using  LGL  points  . 

0/  O/ 

/o - /o 

%Interpolate  a  known  function  f(x)  using  legendre -gauss  -  labatto  points, 
clear 


%Input  Nth  order  interpolation  (N)  and  number  of  evaluation  points  (k) 
N  =  40; 
k  =  50; 

Reinitialization 
errLlnum  =  zeros  (N,  1 )  ; 
errLlden  =  zeros(N,l); 
errLl  =  zeros  (N, 1 )  ; 
errL2num  =  zeros  (N,  1 )  ; 
errL2den  =  zeros  (N, 1); 
errL2  =  zeros (N,l); 
errinfabsn  =  zeros(k,l); 
errinfabsd  =  zeros(k,l); 
errinf  =  zeros(N,l); 

for  n  =  1  :N 

%Setting  up  Grids 
x  =  legendre_gauss_lobatto  (n  +  1) ; 
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29 

31 

33 

35 

37 

39 

41 

43 

45 

47 

49 

51 

53 

55 

57 

59 

61 

63 

65 

67 


z  =  linspace  ( -1  ,1  ,50)  ; 

%Set  up  data  for  Lagrange  Li(Xk) 
f  =  exp( -4*x.A2) ; 

%Constuct  Lagrange 
[L,dL]  =  lagrange_basis (x , z ) ; 

%Evaluation 
Pn  =  f  4-  L ; 

%Error  Analysis 
fex  =  exp ( -4*z . A2)  ; 
for  i  =  l:k 

errLlnum(n)  =  abs  (Pn(  i ) -fex  ( i )  )  +  errLlnum(n)  ; 
errLlden(n)  =  abs ( f ex ( i ) )  +  errLlden(n); 
errL2num(n)  =  (Pn(  i  )  -  fex  ( i )  )  A2  +  errL2num(n)  ; 
errL2den(n)  =  (fex(i))A2  +  errL2den(n)  ; 
errinf absn  ( i )  =  abs (Pn( i ) -fex ( i ) ) ; 
errinfabsd  ( i )  =  abs ( f ex ( i ) )  ; 

end 


errLl(n)  =  errLlnum (n)  / errLl den (n)  ; 
errL2(n)  =  sqrt  (errL2num(n) /errL2den(n)  )  ; 
errinf  (n)  =  max(  errinfabsn  ) /max(  errinfabsd  )  ; 

end 


o/  o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%  1  2  3  4  5  6  7  8 

c  =  [1  0  0;0  0  0;0  0  0;  0  1  0;  0  0  0;  0  0  0;  0  0  0;  1  0  0;... 

0  0  0;0  0  0;0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0.5  0.3  0.2]; 
%  9  10  11  12  13  14  15  16 


0/  o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%plot  Pn  vs.  Actual  Function 

hold  on 

plot (x, f ,  '  *  ) 

t  i  1 1 e  (  ' Lagendre  Gauss  Lobatto  Approximation') 
xlabel (  ' x  ' ) 
ylabel ( ' f (x)  ' ) 
axis ([ -1  1  -0.1  1]) 

legend  (  'N  =  2  '  ,  'N  =  4 ','N  =  8  '  ,  'N  =  16  '  ,  'Exact') 


%Plot  Error  Norms 
figure 
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69 


71 


73 


N  =  1:40; 

semilogy  (N,  errLl  ,N,  errL2  ,N,  errinf) 

t i t 1 e (' Legendre  Gauss  Lobatto  Interpolation  Error') 

xlabel  (  N' ) 

ylabel  ('  Error  Norm') 

legend  (  '  LI  Norm  '  ,  '  L2  Norm  '  ,  '  Inf  Norm  ' ) 


A.2  Integration 


l 


3 


5 


7 


9 


11 


13 


15 


17 


19 


21 


23 


25 


27 


o/  o/ 

/o - /o 

%This  is  the  main  driver  for  Integration  . 

%Written  by  Benjamin  Davis 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%Synopsis  :  Conducting  Integration  using  LGL  points  . 

O/  0/ 

/o - /o 


%Integration  :  Using  the  interpolation  functions  ,  sampling  the 
%basis  functions  at  LGL  and  LG  integration  points  to  perform  the  Gauss 
%quadrature  of  the  given  equation  from  the  project  set  . 
clear 

%Initialization 
N  =  19; 

for  n  =  1  :N 

%Quadratrue  points  and  weights 
[x,w]  =  legendre_gauss_lobatto  (n+1) ; 

%Given  Equation  for  evaluation 
f  =  exp( -4*x.A2)  ; 

%Constuct  Lagrange  matrix  and  differentiation  matrix 
[L,dL]  =  lagrange_basis (x,x) ; 

%Evaluation 
Pn  =  ( f  *L) *  w' ; 

%Evaluation  of  Error 
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29 


31 


33 


35 


37 


exact  =  ( sqrt ( pi ) /2) *  erf (2)  ; 

errLl(n)  =  abs(Pn  -  exact )/ abs ( exact ) ; 

errL2(n)  =  sqrt((Pn  -  exact )  A2/(  exact )  A2) ; 

end 

%Plot  Error 
N  =  1:19; 
semilogy  (N,  errL2 ) 

t i t 1 e (' Legendre  Gauss  Lobatto  Integration  Error') 
xlabel  (  'N' ) 
ylabel ( ' Error ' ) 
legend ('L2  Norm') 
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APPENDIX  B: 

One-Dimensional  Discontinuous  Galerkin 


The  following  are  the  main  MATLAB  codes  for  the  one-dimensional  problem  from 
Chapter  3. 


2 


4 


6 


10 


12 


14 


16 


18 


20 


22 


24 


26 


28 


30 


%This  is  the  Driver  function  using  a  Discontinuous  Galerkin  method  for 
%approximating  a  2nd  Order  Acoustic  Wave  equation  in  one  dimension 
%with  an  Upwind  Flux . 

o/ 

/o 

%Written  by  Benjamin  Davis  Created :  October  2014 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%Synopsis :  Discontinuous  Galerkin  Method  for  wave  equation  in  second 
%order  form  using  inexact  integration  an  Upwind  flux  .  The  outputs 
%are  currently  four  plots  :  Convergence  rates  for  w  and  P  and  plots 
%of  the  numerical  and  exact  solutions  . 

O/  0/ 

clear 

%Initial  Inputs 
N  =  4;  %Polynomial  Order 
n  =  3; 

dtscale  =  1/4; 
c  =  1; 
z  =  1; 

t_final  =  0.58; 
ngl  =  N+l; 

for  Ne  =  2 .  A ( 1 :4 ) 

fprintf  (  'Number  of  Elements  %4d  with  polynomial  order  %2d\n',Ne,N) 
Np  =  Ne*(  ngl ) ; 

%Interpolation  and  Integration  Points 
[psx,w]  =  legendre_gauss_lobatto  ( ngl )  ; 
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34 

36 

38 

40 

42 

44 

46 

48 

50 
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54 
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58 
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62 

64 

66 

68 


%Constuct  Lagrange  matrix  and  differentiation  matrix 
[L,dL]  =  lagrange_basis  (psx  , psx )  ; 

%Construct  Mass  and  Differentiation  Matrices 
[M,D]  =  mass_difflD  (L,  dL,w)  ; 

%Construct  Global  Mass  and  Differentiation  Matrices 
[  coord  ,  intma  ]  =  create_grid  ( ngl  ,Ne,  psx )  ; 
dx  =  coord (N+l) -coord  ( 1 )  ; 

M  =  M.  *dx/2; 

D  =  M\D; 

Dh  =  D'; 

grad2  =  Dh*M*D; 
one  =  ones  (ngl ,  ngl )  ; 
onesM  =  one*M; 


Ml  =  grad2  +  onesM; 


eO  =  sparse (1 ,1 ,1 ,ngl ,1 ) ; 
en  =  sparse  ( ngl ,  1  / 1  /  ngl ,  1 )  ; 


%Building  the  Initial  Condition 
for  e  =  l:Ne 

for  i  =  1  :N+1 

I  =  intma  (e  ,  i ) ; 
x  =  coord ( I ) ; 

q0(e,i,l)  =  n* pi  *  sin  (n* pi * x) ;  %dp/dt=w  at  t=0 
q0(e,i,2)  =  0;  %P  at  t=0 

end 

end 


ql  =  q0(:  ; 

%Time  Step  Calculation 
dt  =  dtscale  *  ( 1  /Ne)  / (NA2) ; 

Nt  =  round ( t_fina  1  /  dt )  ; 
dt  =  t_final/Nt; 

%Periodic  Boundary  Conditions 
[sidep]  =  sideper(Ne); 
%Facemap 

[fmp]  =  facemap(Ne); 


op/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/ffo/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 
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72 

74 

76 

78 

80 

82 

84 

86 

88 

90 

92 

94 

96 

98 

100 

102 

104 

106 

108 


op/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o'o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%Compute  RK  Time- Integration  Coefficients 

kstages  =  4;  %=1  is  RK1 ;  =2  is  RK2;  =  3  is  RK3;  and  =  4  is  RK4 
if (  kstages  <  4) 

[a0,al,beta]  =  compute_ti_coef  ficients  (  kstages  )  ; 


for  k  =  l:Nt 

for  a  =  1 : kstages 

[wf,dpf]  =  upwindflux(Ne,  ql  ,D,  ngl ,  sidep  )  ; 

[R]  =  RHSDG(en ,  eO  ,Dh,D,M,  ql  ,Ne,Ml,wf  ,fmp,  dpf ,  c )  ; 

qm=a0(a) *q0  +  al(a)*ql  +  dt *  beta ( a ) *R; 
ql  =  qm; 

end 

qO  =  qm; 

end 

elseif ( kstages  ==  4) 


a  =  [0,1/2,1/2,11; 
b  =  [1/6,1/3,1/34/6]; 

R  =  0; 

for  k  =  1 : Nt 
qm  =  qO ; 
for  s  =  1:4 

qs  =  qO+dt *a ( s ) *R; 

[wf,dpf]  =  upwindflux  (Ne,  qs  ,D,  ngl  ,  sidep  )  ; 

[R]  =  RHSDG(en ,  eO  ,Dh,D,M,  qs  ,Ne,Ml, wf  ,fmp,  dpf  ,  c  ) ; 

qm  =  qm  +  dt  *b  (  s  )  *R; 

end 

qO  =  qm; 
end 
end 


%Plotting  Purposes 
for  e  =  l:Ne 


for  i  =  1  :N+1 
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132 
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140 

142 
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I  =  intma  ( e  ,  i ) ; 
x  =  coord ( I ) ; 

qex(e,i,l)  =  n*  pi  *  sin  (n*pi  *x)  *  cos  (rupi  *  t_f  inal  )  ; 
qex(e,i,2)  =  sin  (n* pi *x) *  sin  (n* pi  * t_f  inal  )  ; 

end 

end 

%Plotting  purposes 
m=l; 

for  e  =  l:Ne 

for  i  =  1  :N+1 

w(m)  =  qm(  e  ,  i  ,1 )  ; 
wexact(m)  =  qex(e,i,l); 

P (m)  =  qm( e  ,  i  ,2)  ; 
pexact(m)  =  qex(e,i,2); 

I  =  intma  ( e  ,  i ) ; 
xp (m)  =  coord ( I ) ; 
m=3n+ 1 ; 

end 

end 

%Building  the  Error  Plot  Information 
L2errw  =  0; 

L2errp  =  0; 

for  elm  =  l:Ne 

pnts  =  (elm-1)  *(N+1)  +  (1:(N+1)); 
dw  =  w(  pnts  ) -wexact  ( pnts  )  ; 
dp  =  p(pnts  ) -pexact  (pnts  )  ; 

L2errw  =  dw*M*dw'  +  L2errw ; 

L2errp  =  dp*M*dp'  +  L2errp  ; 
end 


L2err (1 ,z ) 
L2err (2 ,z ) 
L2err (3 ,z ) 
L2err (4 ,z ) 
z  =  z  +  1; 

end 


=  Np; 

=  sqrt  (L2errw)  ; 
=  sqrt ( L2errp  ) ; 
=  Ne; 
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148 

150 

152 

154 

156 

158 

160 

162 

164 

166 

168 

170 

172 

174 

176 

178 

180 

182 

184 

186 


O/ 

/o 

°/ 

/o 

C  = 

o/ 

/o 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


1 

[1  0 
0  0 

9 


0;0  1  0;0  0 
0;0  0  0;0  0 
10  11 


3 

0; 

0; 


4 

0  0  0; 
0  0  0; 
12 


5 

0  0 
0  0 
13 


0; 

0; 


6 

0  0 
0  0 
14 


1;  0  0  0; 
0;  0  0  0; 
15 


0.5 


8 

1  0 
0.3 
16 


0;... 

0.21; 


0/  o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 
/o  /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

»W»W«W  Convergence  PI 
figure  ( 1 ) 

subplot  ( 1  ,2 , 1 )  ;  %hold  on 

loglog(L2err(4  ,:)  , L2err ( 2  ,:)  ,  '  *  -  '  ,  'color  '  ,c  (N, : )  ) 

t i 1 1  e  (' Convergence  Plot  for  w' ) 

xlabel(  'Ne  -  Number  of  Elements') 

ylabel  (  'L2  w' ) 

axis  tight 

legend  (  'N  =  4  1  ,  'N  =  6  ','N  =  8  '  ,  'N  =  16  ' ) 

°&8/8/8/8/8/8/8/8/8/8/8/£  Convergence  Plo  1 0/o/8/o/o/o/8/o/o/8/o/o/o/8/o/8/8/8/o/o/8/o/o/o/8/8/o 

subplot (1 ,2 ,2) ;  %hold  on 

loglog(L2err(4  ,:)  ,L2err(3  ,:)  ,  '  *  -  '  ,  'color  '  ,  c  (N, : )  ) 

title  (  'Convergence  Plot  for  p') 

xlabel(  'Ne  -  Number  of  Elements') 

ylabel ( ' L2  p ' ) 

axis  tight 

legend  (  'N  =  4  '  ,  'N  =  6  '  ,  'N  =  8  ','N  =  16  ' ) 

%  %%%%%%%%%%0/o%%%Plot  Numerical  Solution  P(x,  t  )%%%%%%%%%%% 

figure  (2) 

plot  (xp  ,p,xp ,  p exact  ,  '  *  ' ) 

title  (  'Numerical  Solution  for  p(x,t)') 

xlabel ( ' x ' ) 

ylabel  (  '  t  '  ) 

legend  (  '  Numerical  '  ,  '  Exact  ' ) 

%  %%%%%%%%%%%%% Plot  Numerical  Solution  w(x,  t )%%%%%%%%%%% 

figure  (3) 

plot  (xp,w,xp,wexact ,  '  *  ' ) 

title  (  'Numerical  Solution  for  w(x,t)  ') 

xlabel (  'x  '  ) 

ylabel  (  '  t  ' ) 

legend  (  '  Numerical  '  ,  '  Exact  ' ) 

rate_w  =  ( log  (  L2err  (2 , 1  :  end -1 ) )  -  log  (  L2err  (2  ,2  :  end)  )  )  ./  ... 
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188 


190 


192 


194 


( log  (  L2err  (1 ,2  :  end) )  -  log  (  L2err  (1 ,1  :  end  -1) ) )  ; 
rate_p  =  ( log  (  L2err  (3 ,1  :  end -1) ) -log  (  L2err  (3  ,2  :  end) ) )  ./ 
( log  (  L2err  (1 ,2  :  end) )  -  log  (  L2err  (1 ,1  :  end  -1) ) )  ; 

fprintf(  'Convergence  rate  for  w:  ') 
disp  ( rate_w ) 

f printf  (' Convergence  rate  for  p:  ') 
disp  ( rate_p ) 
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/o 


o/ 

/o 


%Function  for  the  Right  Hand  Side  of  Discretized  equations  for 
%one  dimension  DG  from  Chapter  3. 

%Written  by  Benjamin  Davis  Created :  November  2014 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function  [R]  =  RHSDG(en ,  eO  ,Dh,D,M, qn,Ne,Ml,wf  ,fmp,  dpt ,  c ) 
for  e  =  l:Ne 


R(e  , :  ,2)  =  Ml  *  qn  ( e  , :  ,1) 

R(e  , :  ,2 )  =  R(e  , :  ,2)  '  +  Dh*( en*wf  (fmp(2  , e ) )  -  eO *wf  (fmp(  1 , e ) ) ) ; 

R(e , :  , 2 )  =  R(e , : ,2)  '  -  Dh* ( en»en ' » qn(e , :  , 1 )  '  - eO  *e0 ' *  qn(e , :  , 1 )  ' ) 

/ 

R(e,:  ,1)  =  -cA2*Dh*M*D*qn(e,:  ,2) 

R(e,:,l)  =  R(e  , :  ,  1 )  '  +  c  A2t(en»  dpt  (fmp(2  ,e ) )  -eO  *dpf  (fmp(l  ,  e) ) ) ; 

end 

R1  =  M1\R( :  ,2) 

R1  =  R1  ' ; 

R(:  / 2 )  =  Rl; 

R2  =  M\R( :  ,1) 

R2  =  R2  ' ; 

R(:  ,1)  =  R2; 

end 
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/o - /o 

%Function  for  the  Center  Flux  for  the  one  -  dimensional  DG  problem  from 
%Chapter  3 

%Written  by  Benjamin  Davis  Created :  October  2014 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function  [wf,dpf]  =  centerf  lux  (Ne,  qn,D,  ngl  ,  sidep  ) 
for  s  =  l:Ne 

Ls  =  sidep  ( 1 ,  s  )  ; 

Rs  =  sidep (2  ,  s  )  ; 

qLw  =  qn(Ls , ngl ,1) ; 
qRw  =  qn ( Rs  , 1 , 1 ) ; 


dpL  =  D*qn(Ls  , :  ,  2  ) 
dpL  =  dpL(ngl)  ; 
dpR  =  D*  qn  ( Rs  , :  ,  2  )  ' ; 
dpR  =  dpR ( 1 ) ; 

wf  ( s  , : )  =  ( 1 /2 ) * (qLw  +  qRw)  ; 
dpf(s,:)  =  (l/2)*(dpL  +  dpR) ; 

end 

end 
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%Upwind  Flux  function  for  one  dimension  DG.  Used  in  one  dimension  driver 
%for  2nd  Order  Acoustic  wave  equation. 

%Written  by  Benjamin  Davis  Created :  October  2014 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - - - % 


function  [wf,dpf]  =  upwindflux  (Ne,  qn,D,  ngl  ,  sidep  ) 
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c  =  1; 

for  s  =  l:Ne 

Ls  =  sidep  (1 ,  s  )  ; 

Rs  =  sidep  (2  ,  s  )  ; 

qLw  =  qn  ( Ls  ,  ngl  ,  1 )  ; 
qRw  =  qn ( Rs , 1  , 1 )  ; 

dpL  =  D*qn(Ls  , :  ,2 ) 
dpL  =  dpL(ngl); 
dpR  =  D*  qn  ( Rs  , :  ,  2 )  ' ; 
dpR  =  dpR(l)  ; 

wf(s,:)  =  (1/2)* (qLw  +  qRw)  +  c/2*(dpR  -  dpL)  ; 
dpf  ( s  , : )  =  (l/2)*(  dpL  +  dpR )  +  1  /  ( 2  *  c )  *  (qRw  -  qLw) ; 

end 

end 
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0/  o/ 

/o - /o 

%Code  given  to  Professor  's  F.X.  Giraldo  MA4245  Class  July  2014 
%Used  by  Ben  Davis 

%This  code  computes  the  Legendre -Gauss - Lobatto  points  and  weights 
%which  are  the  roots  of  the  Lobatto  Polynomials . 

%Written  by  F.X.  Giraldo  on  4/2000 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 


function  [xgl,wgl]  =  legendre_gauss_lobatto  (P) 
p=P - 1 ; 

ph=floor  (  (p  +  l)/2  ); 
for  i=l:ph 

x=cos(  (2*  i -1 )  *  pi /( 2*p  +  l)  ); 
for  k  =  l:20 

[L0,L0_1  ,L0_2]  =  legendre_poly  (p,x) ; 
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%Get  new  Newton  Iteration 

dx  =  -(l-xA2)  *L0_1  /(  -2 * x*L0_l  +  (1  -xA2)  *L0_2 ) ; 
x=x+dx ; 

if  (abs(dx)  <  1.0e-20) 
break 

end 

end 

X§1 (p+2-i )=x; 

wgl(p+2-i )  =  2/(p*(p  +  l)  *L0 A2) ; 

end 

%Check  for  Zero  Root 
if  (p+1  ~=  2*ph) 
x  =  0; 

[L0,L0_1  ,L0_2]  =  legendre_poly  (p,x) ; 
xgl  (ph  +  l)=x; 

wgl(ph  +  l)  =  2/(p*(p  +  l)*L0A2) ; 

end 

%Find  remainder  of  roots  via  symmetry 
for  i=l:ph 

xgl(i)=-xgl(p+2-i); 
wgl(  i  )=+wgl(p+2-i )  ; 

end 
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%Function  for  building  the  Lagrange  Polynomials . 

%Written  by  Benjamin  Davis  in  MA4245  Created :  July  2014  in  MA4245 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function  [L,dL]  =  lagrange_basis  (x,z) 

%Nth  order  interpolation 
n  =  length (x) ; 

%Length  of  the  equally  spaced  grid  for  k  =  1:50 
h  =  length ( z )  ; 
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%Initialize  the  Lagrange  Matrix 
L  =  ones (n,h) ; 
dL  =  zeros (n,h) ; 

%Computation  for  Lagrange  Matrix 
for  k  =  1 : h 

for  i  =  l:n 

for  j  =  l:n 
dl  =  1; 

if  j  ~=  i  %  If  j  does  not  equal  i 

%Equation  for  the  Lagrange  Polynomial 
L(i,k)  =  (z(k)-x(j)) ./(x(i)-x(j))  *  L(i,k); 
for  1  =  l:n 

if  (1  ~=  i)  &&  (1  ~=  j ) 

dl  =  dl*(z(k)-x(l)) ./(x(i)-x(l)); 

end 

end 

dL(i,k)  =  dL(i,k)  +  dl /  ( x ( i ) -x( j ) ) ; 

end 

end 

end 

end 

end 
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/o - /o 

%Code  given  by  Professor  F.X.  Giraldo  to  MA4245  class. 

%Used  by  Ben  Davis 

%This  code  computes  the  Legendre  Polynomials  and  its  1st  and  2nd 
%d  erivatives 

%Written  by  F.X.  Giraldo  on  4/2000 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%This  code  was  written  by  Professor  Giraldo  and  given  to  his  MA4245 
%Galerkin  Methods  class  in  July  2014. 

O/  O/ 

/o - /o 
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function  [L0 ,  L0_1 ,  L0_2  ]  =  legendre_poly  (p,x) 

L1=0;L1_1=0;L1_2=0; 

LO  =  1;L0_1  =  0;L0_2  =0; 

for  i=l:p 

L2=L1 ; L2_1=L1_1 ; L2_2=L1_2 ; 

L1=L0 ; L1_1=L0_1 ; L1_2=L0_2 ; 
a  =(2*  i  -1)  /  i ; 

b=(i  -1)  /  i ; 

L0=a*x*Ll  -  b*L2; 

L0_l=a*(Ll  +  x*Ll_l )  -  b*L2_l ; 
L0_2=a*(2*Ll_l  +  x*Ll_2)  -  b*L2_2; 

end 
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%Function  for  building  the  Mass  and  Differentiation  matrices  for 
%Qne- Dimension .  Used  in  Thesis  ID  Upwind  code. 

%Written  by  Benjamin  Davis  Created:  July  2014 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  % 

function  [M,D]  =  mass_difflD  (L,dL,w) 
n  =  length (L ( : ,1) ) ; 

M  =  zeros (n,n) ; 

D  =  zeros (n,n) ; 

for  i  =  l:n 


for  j  =  l:n 

M(i,j)  =  ((L(i  ,:)  .*L(  j  , : )  )  *w  ( 1  ,:)  ')  ; 
D(i,j)  =  ((L(i  ,:)  ■  * dL ( j  ,:)  ).w(l  ,:)  '); 

end 

end 

end 
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%Code  given  to  Professor  's  F.X.  Giraldo  's  MA4245  class 
%Used  by  Ben  Davis 

%This  function  computes  the  LGL  grid  and  elements . 
%Written  by  F.X.  Giraldo  on  10/2003 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 


function  [ coord  ,  intma ]  =  create_grid  (ngl  ,nelem  ,  xgl ) 
%Set  some  constants 
xmin  =  - 1 ; 
xmax=  +  l; 

dx  =  (xmax-xmin)  /nelem ; 

%Generate  Grid  Points 

ip  =1; 

coord  ( 1 )  =xmin ; 
for  i=l:nelem 

x0=xmin  +  (i-l)*dx; 
intma  ( i  ,  1 )  =ip  ; 
for  j  =2:  ngl 
ip=ip  +  1; 

coord(ip)=(  xgl(j)+l  )*dx/2  +  xO ; 
intma  ( i  ,  j  )=ip  ; 

end 

end 


2  %  Function  for  implementing  periodic  boundary  conditions  for  Thesis  ID 
%  Upwind  code  . 

4  %Written  by  Benjamin  Davis  Created :  November  2014 
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0/ 

/o 

Department  of  Applied  Mathematics 

0/ 

/o 

Naval  Postgraduate  School 

O/ 

/o 

Monterey,  CA  93943-5216 

0/ 

/o - 

O/ 

- /o 

function  [sidep]  =  sideper(Ne) 
sidep  =  zeros(2,Ne); 
for  e  =  l:Ne 

sidep  (1  ,e)  =  e -1; 
sidep  (2  ,e)  =  e; 
if  e==l 

sidep  ( 1 ,  e )  =  Ne; 

end 

end 

end 
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%Facemaping  periodic  function  for  1 -Dimension  and  used  in  Thesis  ID 
%Upwind  code . 

%Written  by  Benjamin  Davis  Created :  November  2014 


O/ 

/o 

Department  of  Applied  Mathematics 

0/ 

/o 

Naval  Postgraduate  School 

0/ 

/o 

Monterey,  CA  93943-5216 

0/ 

/o - 

O/ 

- /o 

function  [fmp]  =  facemap(Ne) 
fmp  =  zeros  (2, Ne); 
for  e  =  l:Ne 

fmp  ( 1  ,  e )  =  e ; 
fmp  ( 2  ,  e )  =  e  + 1 ; 
i  f  e==Ne 

fmp  ( 2  ,  e )  =  fmp  (1,1); 

end 

end 

end 
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APPENDIX  C: 

Two-Dimensional  Discontinuous  Galerkin 


The  following  are  the  main  MATLAB  codes  for  the  two-dimensional  problem  from 
Chapter  4. 
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0/  o/ 

%This  is  the  main  driver  for  approximating  the  2D  Acoustic  Wave 
%Equation  with  an  upwind  flux  on  a  washer  grid  .  The  washer  grid 
%can  be  scaled  to  various  sizes  based  off  of  user  input. 

%Written  by  Benjamin  Davis  and  Asst.  Professor  Jeremy  Kozdon 
%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%Synopsis :  Discontinuous  Galerkin  Method  for  wave  equation  in  second 
%order  form  using  inexact  integration  and  an  upwind  flux  .  The  outputs 
%are  currently  five  plots  :  Convergence  rates  for  w  and  P  and  plots  of 
%the  numerical  and  exact  solutions  . 

O/  O/ 

/o - /o 

clear 


N  =  2; 
n  =  1; 
c  =  1; 

t_final  =  pi/2; 

%Define  beta  to  be  zero  and  will  run  with  no  theta  dependence, 
beta  =  4; 

%  skew_mesh  =  0  :  :  rectangular  mesh 

%  skew_mesh  =1  ::  skew  element  (straight  sided) 

%  skew_mesh  =  2  ::  skew  element  (curved  elements) 

skew_mesh  =  2 ; 

%Define  anonymous  functions  for  the  solution 

r_ex  =  @(x,y)  sqrt(xA2  +  yA2) ; 

theta  =  @(x,y)  acos (x/r_ex (x,y) ) ; 

g_ex  =  @(r_ex)  besselj (beta  ,r_ex) ; 

f_ex  =@(x,y,t)  sin  ( c *t - beta  *  theta (x ,y) ) ; 
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P_ex  =@(x,y,t)  f_ex(x,y, t ) *g_ex( r_ex (x,y) ) ;  %p(x,y,t) 

w_ex  =@(x,y,t)  c  *  cos  (  c  *t  -  beta  *  theta  (x  ,y)  )  *g_ex  (  r_ex  (x  ,y)  )  ;  %dp/dt 

%Select  inner  and  outer  radius  of  washer. 

%Found  these  values  using  fzero  command  on  Bessel  function . 
rl  =  fzero  ( g_ex  ,6 )  ; 
r2  =  fzero  ( g_ex  ,  1 4)  ; 
disp  ([  rl  ,  r2  ])  ; 

%Full  Washer. 

%rn  is  what  2*pi  will  be  divided  by  to  change  the  size  of  the  mesh. 
%For  example,  if  you  use  4,  you  will  get  pi/2  or  1/4  of  the  washer, 
rn  =  4; 

%Scaling  the  washer. 

%rn  =  ceil  ((2*pi*r2) /(  r2 -rl  ) )  ; 

%For  storing  information  ,  Don' t  Change 
z  =  l; 

for  nel  =  2.A(1:3) 

nq  =  N+l;  “/Number  of  integration  /  quadrature  points 

ngl  =  N+l;  “/Number  of  interpolation  points  in  one  direction 

nelx  =  nel ; 

nely  =  nel ; 

“/.Interpolation  and  Integration  Points 
[psx,w]  =  legendre_gauss_lobatto  (N+l) ; 

“/.Create  Grid 

[  coord  ,  intma  ,bsido  ,  iperiodic  ,Np,Ne,  nboun  ,nface]  =  create_grid_2d( 
nelx  ,  nely  ,N,  psx  ,  plot_grid  ,  skew_mesh) ; 

fprintf  (  'Number  of  Elements  %4d  with  polynomial  order  %2d\n',Ne,N) 

%Create  Sides /Edge  Information  for  EG 

[iside,jeside]  =  create_side(  intma  ,  bsido  ,  Np,Ne,  nboun,  nf  ace  ,  ngl )  ; 

[  face  ,  imapl ,  imapr  ]  =  create_f  ace  ( iside  ,  intma  ,  nface  ,  ngl )  ; 
s_face  =  face ; 

“/Changing  face  code  to  enforce  new  BC. 

[face]  =  faceBC  ( nface  ,  face  )  ; 
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%Will  not  be  periodic  with  face  from  faceBC .  If  you  want  periodic, 
%comment  out  faceBC  function . 

face  =  create_face_periodicity  ( iside  ,  face  ,  coord  ,  nface  ,nboun)  ; 
%Building  the  washer. 

[Rad]  =  radius ( coord  ,  rl  ,  r2  )  ; 

[theta]  =  thetapolar ( coord , rn )  ; 

[coord]  =  newcoord  (Rad,  theta  ,  coord  )  ; 

%Construct  Lagrange  Basis  and  Jacobian  Matrix 
[L,dL]  =  lagrange_basis  (psx  ,  psx )  ; 

[  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  ,  x_ksi  ,  x_eta  ,  y_ksi  ,  y_eta  ,  jac  ]  =  metrics2  ( 
coord  ,  intma  ,L,dL,Ne,  ngl  ,nq)  ; 

%Building  and  Element  to  Element  function  . 

%EtoEfunctBC  incorporates  boundary  conditions  . 

EtoEBC  =  EtoEfunctBC  ( nface  ,  face  ,Ne) ; 

%Constuct  M,  D,  D_ksi ,  D_eta  Matrices 
[M,D]  =  mass_diff2D  (L,  dL,w)  ; 

M_1D  =  M; 

M  =  kron  (M,M) ; 

D_eta  =  kron  (D,  eye  ( ngl )) ; 

D_ksi  =  kron  ( eye  ( ngl )  ,D) ; 

%Construct  L1,L2,L3,L4 
eO  =  sparse (1 ,1 ,1 ,ngl ,1 )  ; 
en  =  sparse(ngl,l ,1 ,ngl,l) ; 

LI  =  kron ( eO ', eye (ngl )) ; 

L2  =  kron  ( eye  ( ngl )  ,en ')  ; 

L3  =  kron ( en ', eye (ngl )) ; 

L4  =  kron ( eye ( ngl ) , eO ' ) ; 

%Building  Matrix  Terms  and  new  Jacobian  for  solving  the  RHS  for 
equation  5,  6  and  7. 

[MAT, MAT_ksi , MAT_eta ,  J e  ,A,B,C,D,E,F,G,H]  =  MatrixTerms2D (Ne,  D_ksi , 
D_eta  ,  y_eta  ,  y_ksi  ,  x_eta  ,  x_ksi  ,  jac  ,M) ; 

%Building  Dx  and  Dy  for  RHS 

[Dx,Dy]  =  DxDy(Ne,  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  ,  D_ksi ,  D_eta ) ; 

%Building  the  surface  jacobians  for  the  faces 

[  Sjl  ,  Sj 2  ,  Sj 3  ,  Sj 4  ]  =  SurfaceJac2D  (Ne,  LI  ,L2  ,L3 ,  L4 ,  x_ksi  ,  x_eta  ,  y_ksi  , 
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y_eta ) ; 

%Compute  the  Normals  per  element 

[nx_l  ,  nx_2  ,  nx_3  ,  nx_4  ,  ny_l  ,  ny_2  ,  ny_3  ,  ny_4  ]  =  ElementNormals2D  (Ne,  LI , 
L2  ,L3 ,  L4 ,  x_ksi  ,  x_eta  ,  y_ksi  ,  y_eta  ,  Sjl  ,  Sj2  ,  S j  3  ,  S j  4  ) ; 

%Big  Matrix  of  Ones 
Ml  =  ones  ( ngl  *  ngl )  ; 

%Building  Initial  Condition 
for  e  = 
for 


1  :Ne 

i  =  1: 

ngl 

for  j 

=  1 :  ngl 

I 

=  intma(e 

,i ,  j ); 

X 

=  coord (I 

,1); 

y 

=  coord (I 

,2); 

q0 

{ 1  }(e,  i  ,  j 

)  =  w_ex  ( x ,  y ,  0 ) 

q0{2}(e,i,j)  =  P_ex(x,y,0) 


%dp/dt  =  w  at  t=0 

%P(x,y,t)  at  t=0; 


end 

end 

end 


ql  =  qO; 


%  Time  Step  Calculation 
dt  =  (1/2)  *((1/Ne)/(NA2)); 
Nt  =  round ( t_fina  1  /  dt )  ; 
dt  =  t_final/Nt; 


%  Beginning  of  the  RK54  Time  integrator 
a  =  [  0.0, 

-567301805773.0/  1357537059087.0 , 
-2404267990393.0/  2016746695238.0, 
-3550918686646.0/  2091501179385.0, 
-1275806237668.0/  842570457699.0]; 


b  =  [1432997174477.0/ 
5161836677717.0/ 
1720146321549.0/ 
3134564353537.0/ 


9575080441755.0, 

13612068292357.0, 

2090206949498.0, 

4481467310338.0, 
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176 


2277821191437.0/  14882151754819.0]; 


R{ 1 1  =  0; 

R{2}  =  0; 

for  k  =  l:Nt 
qm  =  qO ; 
for  s  =  1:5 

R{ 1 }  =  a( s ) *R{  1 ) ; 

R{ 2 }  =  a(s)*R{2}; 

%Building  all  the  flux  information  with  Boundary  Conditions 
[fwl ,  fw2 ,  fw3  , fw4 ,  fpl  ,  fp2  ,  fp3  ,  fp4  ]  =  UpwindFlux2DFwFpBC  (Ne,qm 
,  LI , L2 , L3 , L4 ,  EtoEBC  ,Dx,Dy,  nx_l  ,  ny_l ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4 )  ; 
%Solving  Right  Eland  Side 

[Rl]  =  RFlS]Xi2]>i(Ne,MAT,Je  ,M,M_lD,Ml,qm,Dx,Dy,Ll,L2,L3,L4, 
nx_l  ,  ny_l  ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4  ,  Sj  1  ,  Sj  2  ,  Sj  3  ,  Sj 4  ,  fwl ,  fw2  ,  fw3  , 
fw4 ,  fp  1  ,  fp2  ,  fp3  ,  fp4  ,  c  ,  MAT_ksi ,  MAT_eta ) ; 

R{ 1 }  =  Rl { 1 }  +  R{ 1 } ; 

R{ 2 }  =  Rl { 2 }  +  R{ 2  } ; 

qm{l}  =  qmjl}  +  dt*b(  s  )  *R{  1 } ; 
qm{2}  =  qm{2}  +  dt*b(  s  )  *R{  2  } ; 

end 

qO  =  qm; 

end 


%For  plotting  purposes 

%Solving  the  exact  solution  at  final  time, 
for  e  =  l:Ne 

for  i  =  1 : ngl 

for  j  =  1 : ngl 
I  =  intma  (e  ,  i  ,  j  )  ; 
x  =  coord (1,1) ; 
y  =  coord (1,2) ; 

qex  { 1 }  (  e  ,  i  ,  j  )  =  w_ex  (x,y,  t_final  )  ; 
qex { 2  }  ( e , i , j )  =  P_ex(x,y, t_final ) ; 
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end 


end 

end 

%Plotting  the  exact  solution  vs.  numerical  solution. 
m=l; 

for  e  =  l:Ne 

for  i  =  l:ngl 

for  j  =  1 :  ngl 

w(m)  =  qm{  1 }  (e,  i  ,  j  )  ; 
wexact(m)  =  qex  { 1 }  ( e ,  i  ,  j  )  ; 

p(m)  =  qm{2}(e,i  ,  j  ); 
pexact(m)  =  qex { 2  }  ( e ,  i  ,  j  )  ; 

I  =  intma(e,  i  ,  j  )  ; 
xp(m)  =  coord ( I ) ; 
m=m+l; 

end 

end 

end 

%Building  the  Error  Plot  Information 
L2errw  =  0; 

L2errp  =  0; 

for  elm  =  l:Ne 

pnts  =  (elm-1)  *(N+1)A2  +  (1:(N+1)A2); 
dw  =  w(  pnts  )  -wexact  (  pnts  )  ; 
dp  =  p(  pnts  ) -pexact  ( pnts  )  ; 

L2errw  =  dw*M*  Je  { elm }  *dw'  +  L2errw ; 

L2errp  =  dp*M*  Je  { elm }  *  dp'  +  L2errp  ; 

end 

L2err(l,z)  =  sqrt(Np); 

L2err(2,z)  =  sqrt  (L2errw)  ; 

L2err(3,z)  =  sqrt ( L2errp  ) ; 
z  =  z  +  1; 

%Convergence  Rates 

rate_w  =  ( log  (  L2err  (2 , 1  :  end -1 )  ) -log  (  L2err  (2  ,2  :  end)  )  )  ./  (log(L2err 
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(1,2:  end)  )-log(L2err(l,l  :  end  -1 ) )  )  ; 

rate_p  =  ( log  (  L2err  (3 , 1  :  end -1 )  ) -log  (  L2err  (3  ,2  :  end)  )  )  ./ 
(1,2:  end)  )-log(L2err(l,l  :  end  -1 ) )  )  ; 

fprintf  (' Convergence  rate  for  w:  ') 
disp  ( rate_w ) 

fprintf (' Convergence  rate  for  p:  ') 
disp  ( rate_p  ) 

end 


oa/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/cyo/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

O/J/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/ 
Ay o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


o/ 

/o 

c  = 


1 

[1  0 
0  0 

9 


0;0  1 
0;0  0 
10 


0;0  0 
0;0  0 
11 


3  4  5  6  7  8 

0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0.1  0.2 

0;  0  0  0;  0  0  0;  0  0  0;  0.5  0.3 

13  14  15  16 


0;  0  0 

%  9  10  11  12 

figure  (1)  ; 

w««mw  Convergence  Pi  O  fflo/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

subplot  ( 1  ,2 , 1 )  ;  %hold  on 

loglog(L2err(l  ,:)  ,L2err(2  ,:)  ,  '  *  -  '  ,  'color  '  ,  c  (N, : )  ) 

t i 1 1  e  (' Convergence  Plot  for  w' ) 

xlabel (  'Np  ' ) 

ylabel  (  'L2  w' ) 

axis  tight 

»«««%V  Convergence  Plo  ^C/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

subplot (1 ,2 ,2)  ;%  hold  on 

loglog(L2err  (1  ,:)  ,L2err(3  ,:)  ,  '  *  -  '  ,  'color  '  ,  c  (N, : )  ) 

t i 1 1  e  (' Convergence  Plot  for  p') 

xlabel (  'Np  ' ) 

ylabel ( ' L2  p ' ) 

axis  tight 

%Plot  Numerical  Solution  P(x,y,t) 
figure  (2)  ; 

t  i  1 1  e  (' Numerical  Solution  for  P') 
xlabel (  x  ' ) 
ylabel  (  ' y  ' ) 
zlabel (  '  t  ' ) 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  ,  2 )  ,  qm{  2  }(:),'*’ ) 


( log  ( L2err 


0.8;... 

0.2]; 
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%Plot  Exact  Solution  of  P(x,y,t) 
figure  (3)  ; 

title  (  Exact  Solution  for  P') 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  , 2 )  ,  qex  {  2  }(:),'*' ) 
%Plot  Numerical  Solution  w 
figure  (4)  ; 

t  i  1 1  e  (  ' Numerical  Solution  for  w' ) 

%  xlabel ( ' x ' ) 

%  ylabel ( 'y ' ) 

%  z lab el ( ' t  ' ) 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  , 2 )  ,  qm{  1  }(:),'*' ) 
%Plot  Exact  Solution  of  w 
figure  (5)  ; 

title  (  Exact  Solution  for  w') 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  , 2 )  ,  qex  {  1  }(:),'*' ) 
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%This  is  the  main  driver  for  approximating  the  2D  Acoustic  Wave 
%Equation  with  Upwind  Flux  on  a  Square  Grid.  The  grid  can  be  rotated 
%various  degrees  and  skewed  based  on  user  input . 

%Written  by  Benjamin  Davis  and  Asst.  Professor  Jeremy  Kozdon 
%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

°/<%Synopsis :  Discontinuous  Galerkin  Method  for  wave  equation  in 

%second  order  form  using  inexact  integration  and  upwind  flux  .  The 
%outputs  are  currently  four  plots  :  Convergence  rates  for  w  and  P  and 
%plots  of  the  numerical  and  exact  solutions  . 

O/  0/ 

/o - /o 

clear 

N  =  6; 
n  =  1; 
c  =  1; 

t_final  =  .25; 
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%  skew_mesh  =  0  : :  rectangular  mesh 

%  skew_mesh  =1  ::  skew  element  (straight  sided) 

%  skew_mesh  =  2  ::  skew  element  (curved  elements) 

skew_mesh  =  2 ; 

%Rotation  Angle 

grid_rotation_angle  =0;  °/£3CW  rotation  in  degrees 
%For  information  storage.  Don't  change, 
z  =  1; 

for  nel  =  2.A(1:3) 

nq  =  N+l;  %Number  of  integration  /  quadrature  points 
ngl  =  N+l;  %Number  of  interpolation  points  in  one  direction  of  an 
element 
nelx  =  nel ; 
nely  =  nel ; 

plot_grid  =  0;  %1  will  display  the  grid,  0  will  not  display  the  grid 

%Interpolation  and  Integration  Points 
[psx,w]  =  legendre_gauss_lobatto  (N+l) ; 

%Create  Grid 

[ coord  , intma  ,bsido  ,  iperiodic  ,Np,Ne, nboun ,nface]  =  create_grid_2d( 
nelx  ,  nely  ,N,  psx  ,  plot_grid  ,  skew_mesh) ; 

fprintf  (  'Number  of  Elements  %4d  with  polynomial  order  %2d\n',Ne,N) 
%Rotate  Grid 

[coord_rotated]  =  rotate_grid_v2  ( coord  ,  intma  ,Np,Ne,  ngl  ,plot_grid  , 
grid_rotation_angle  )  ; 

%Store  Rotated  COORDS 
coord=coord_rotated  ; 

%Construct  Lagrange  Basis  and  Jacobian  Matrix 
[L,dL]  =  lagrange_basis  (psx  ,  psx )  ; 

[  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  ,  x_ksi  ,  x_eta  ,  y_ksi  ,  y_eta  ,  jac  ]  =  metrics2  ( 
coord  ,  intma  ,L,dL,Ne,  ngl  ,nq)  ; 

%Create  Sides /Edge  Information  for  EG 

[iside,jeside]  =  create_side(  intma  ,  bsido  ,  Np,Ne,  nboun,  nf  ace  ,  ngl )  ; 

[  face  ,  imapl ,  imapr  ]  =  create_f  ace  ( iside  ,  intma  ,  nface  ,  ngl )  ; 
face  =  create_face_periodicity  ( iside  ,  face  ,  coord  ,  nface  , nboun)  ; 
%Building  and  Element  to  Element  function 
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EtoE  =  EtoEfunct  ( nface  ,  face  ,Ne)  ; 

%Constuct  M,  D,  D_ksi ,  D_eta  Matrices 
[M,D]  =  mass_diff2D  (L,  dL,w)  ; 

M_1D  =  M; 

M  =  kron  (M,M)  ; 

D_eta  =  kron  (D,  eye  ( ngl )) ; 

D_ksi  =  kron  ( eye  ( ngl )  ,D) ; 

%Construct  L1,L2,L3,L4 
eO  =  sparse  (1  ,1  ,1  ,ngl  ,1 )  ; 
en  =  sparse(ngl,l  ,l,ngl,l)  ; 

LI  =  kron  ( eO ',  eye  (ngl )) ; 

L2  =  kron ( eye ( ngl ) ,en ') ; 

L3  =  kron ( en eye (ngl )) ; 

L4  =  kron ( eye ( ngl ) , eO ' ) ; 

%Building  Matrix  Terms  and  new  Jacobian  for  solving  the  RHS  for 
equation  5,  6  and  7. 

[MAT, MAT_ksi , MAT_eta ,  Je  , A, B,C,D,E, F ,G,H]  =  MatrixTerms2D (Ne,  D_ksi , 
D_eta  ,  y_eta  ,  y_ksi  ,  x_eta  ,  x_ksi  ,  jac  ,M) ; 

%Building  Dx  and  Dy  for  RHS 

[Dx,Dy]  =  DxDy(Ne,  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  ,  D_ksi ,  D_eta ) ; 

%Building  the  surface  jacobians  for  the  faces 

[  Sjl  ,  S j 2  ,  Sj 3  ,  S j 4  ]  =  SurfaceJac2D  (Ne,  LI  ,L2  ,L3 ,  L4 ,  x_ksi  ,  x_eta  ,  y_ksi  , 
y_eta ) ; 

%Compute  the  Normals  per  element 

[nx_l  ,  nx_2  ,  nx_3  ,  nx_4  ,  ny_l  ,  ny_2  ,  ny_3  ,  ny_4  ]  =  ElementNormals2D  (Ne,  LI , 
L2  ,L3 , L4 ,  x_ksi  ,  x_eta  ,  y_ksi  ,  y_eta  ,  Sjl  ,  Sj 2  ,  Sj 3  ,  Sj 4  ) ; 

%Big  Matrix  of  Ones 
Ml  =  ones  ( ngl  *  ngl )  ; 

%Building  Initial  Condition 

%Equation  for  initial  condition  is  the  following: 

%dp/dt  =  n*pi  *  sin  (n*  pi  *x)  *  sin  (n*pi  *y )  *  cos  (n*  pi  *  t ) 

%P(x,y, t)  =  sin(n*pi*x)*sin(n*pi*y)*sin(n*pi*t) 
for  e  =  l:Ne 

for  i  =  1 : ngl 

for  j  =  1 : ngl 

I  =  intma  (e  ,  i  ,  j  )  ; 
x  =  coord (1,1) ; 
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y  =  coord  (1,2)  ; 

q0{l)(e,i,j)  =  sqrt (2 ) *n* pi  *  sin  (n* pi *x) *  sin  (n* pi *y ) ; 
q0 { 2 } ( e  ,  i  ,  j  )  =  0; 

end 

end 

end 


ql  =  qO; 


%  Time  Step  Calculation 
dt  =  ( ( 1  /Ne)  / (NA2) )  ; 

Nt  =  round ( t_fina  1  /  dt )  ; 
dt  =  t_final/Nt; 


%  Beginning  of  the  RK54  Time  integrator 
a  =  [  0.0, 

-567301805773.0/  1357537059087.0 , 

-2404267990393.0/  2016746695238.0, 

-3550918686646.0/  2091501179385.0, 

-1275806237668.0/  842570457699.0]; 

b  =  [1432997174477.0/  9575080441755.0, 

5161836677717.0/  13612068292357.0, 

1720146321549.0/  2090206949498.0, 

3134564353537.0/  4481467310338.0, 

2277821191437.0/  14882151754819.0]; 

R{ 1 }  =  0; 

R{2]  =  0; 

for  k  =  1 : Nt 
qm  =  qO ; 
for  s  =  1:5 

R{ 1 }  =  a( s ) *R{  1 } ; 

R{ 2 }  =  a( s ) *R{ 2 } ; 

%Building  all  the  flux  information 

[fwl ,  fw2 ,  fw3  , fw4 ,  fpl  ,  fp2  ,  fp3  ,  fp4  ]  =  UpwindFlux2DFwFp(Ne,qm, 
LI ,  L2  ,  L3  ,  L4 ,  EtoE  ,Dx,Dy,  nx_l  ,  ny_l  ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4 )  ; 
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%Solving  Right  hand  side 

[Rl]  =  RHSDG2Dn ( Ne ,MAT,  J e  ,M, M_1D, Ml, qm, Dx , Dy , LI , L2 , L3 , L4 , 
nx_l ,  ny_l ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4  ,  Sj  1  ,  Sj  2  ,  S j  3  ,  Sj  4  ,  fwl ,  fw2 ,  fw3 


,  fp2  , 

fp3  ,  fp4  , 

c  ,  MAT_ksi ,  MAT_eta ) 

R{  l } 

=  Rl  { 1 } 

+  R{  1 }; 

R{  2 } 

=  Rl  { 2 } 

+  R{ 2 ) ; 

qmjl 

}  =  qm{i] 

l  +  dt*b( s ) *R{ 1 } ; 

qm{2 

}  =  qm{2] 

l  +  dt*b( s ) *R{ 2 } ; 

end 

qO  =  qm; 

end 

0/3/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/ 

/ao/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

O/D/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/ 

/o'o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

%For  plotting  purposes 

%Solving  the  exact  solution  at  final  time, 
for  e  =  l:Ne 

for  i  =  1 : ngl 

for  j  =  1 : ngl 
I  =  intma  (e  ,  i  ,  j  )  ; 
x  =  coord (1,1) ; 
y  =  coord (1,2) ; 

qex { 1 } ( e , i , j )  =  sqrt (2 ) *n* pi  *  sin  (n*pi *x ) *  sin  (n* pi *y ) *  cos ( 
sqrt(2)*n*pi*t_final) ; 

qex  {  2  }  (  e  ,  i  ,  j  )  =  sin  (n*  pi  *x)  *  sin  (n*pi  *y )  *  sin  (  sqr  t  ( 2 )  *n*  pi  * 

t_f inal ) ; 

end 

end 

end 

%Plotting  the  exact  solution  vs.  numerical  solution. 
m=  1 ; 

for  e  =  l:Ne 

for  i  =  1 : ngl 

for  j  =  1 : ngl 

w(m)  =  qm{  1 }  (e  ,  i  ,  j  )  ; 
wexact(m)  =  qex  { 1 }  ( e ,  i  ,  j  )  ; 
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P  (m)  =  qm{2}(e,i  ,  j  ); 
pexact(m)  =  qex { 2  }  ( e ,  i  ,  j  )  ; 


I  =  intma(e,  i  ,  j  )  ; 
xp(m)  =  coord ( I ) ; 
rrmn+1; 

end 

end 

end 

%Building  the  Error  Plot  Information 
L2errw  =  0; 

L2errp  =  0; 


for  elm  =  l:Ne 

pnts  =  (elm-1)  *(N+1)A2  +  (1:(N+1)A2); 
dw  =  w(  pnts  ) -wexact  (  pnts  )  ; 
dp  =  p(  pnts  ) -pexact  ( pnts  )  ; 

L2errw  =  dw*M*  Je  { elm  }  *dw'  +  L2errw ; 
L2errp  =  dp*M*  Je  { elm }  *  dp'  +  L2errp  ; 


end 

L2err(l,z)  =  sqrt(Np); 

L2err(2,z)  =  sqrt  (L2errw)  ; 

L2err(3,z)  =  sqrt ( L2errp  ) ; 
z  =  z  +  1; 

%Convergence  Rates 

rate_w  =  ( log  (  L2err  (2 , 1  :  end -1 )  ) -log  (  L2err  (2  ,2  :  end)  )  ) 
(1,2:  end)  )-log(L2err(l,l  :  end  -1 ) )  )  ; 

rate_p  =  ( log  (  L2err  (3 , 1  :  end -1 )  ) -log  (  L2err  (3  ,2  :  end)  )  ) 
(1,2:  end)  )-log(L2err(l,l  :  end  -1 ) )  )  ; 


./  (log(L2err 
./  (log(L2err 


fprintf(  'Convergence  rate  for  w:  ') 
disp  ( rate_w ) 


fprintf  (' Convergence  rate  for  p:  ') 
disp  ( rate_p  ) 

end 


%  1  2  3  4  5  6  7  8 
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c  =  [1  0  0;0  1  0;0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  ( 
0  0  0;0  0  0;0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  ( 
%  9  10  11  12  13  14  15 

figure  (1)  ; 

Convergence  Plof/o%%%%%%%%%%%0/o%%% 

subplot (1 ,2 ,1 ) ;  %hold  on 

loglog(L2err(l  ,:)  ,L2err(2  ,:)  ,  'color  '  ,  c  (N, : )  ) 

t i 1 1  e  (' Convergence  Plot  for  w' ) 

xlabel (  'Np  ' ) 

ylabel  (  'L2  w' ) 

axis  tight 

»«««°»  Convergence  Plof/o%%%%%%%%%%%%%% 

subplot (1 ,2 ,2)  ;%  hold  on 

loglog(L2err(l  ,:)  ,L2err(3  ,:)  ,  'color  '  ,  c  (N, : )  ) 

title  (  'Convergence  Plot  for  p') 

xlabel (  'Np  ' ) 

ylabel ( ' L2  p ' ) 

axis  tight 

%Plot  Numerical  Solution  P(x,y,t) 
figure  (2)  ; 

t  i  1 1  e  (' Numerical  Solution  for  P') 
xlabel (  x  ' ) 
ylabel  (  ' y  ' ) 
zlabel (  '  t  ' ) 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  , 2 )  ,  qm{  2  }(:),'*' ) 
%Plot  Exact  Solution  of  P(x,y,t) 
figure  (3)  ; 

title  (  Exact  Solution  for  P') 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  ,  2 )  ,  qex  {  2 }(:),'*' ) 
%Plot  Numerical  Solution  w 
figure  (4)  ; 

t  i  1 1  e  (' Numerical  Solution  for  w' ) 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  , 2 )  ,  qm{  1  }(:),'*' ) 
%Plot  Exact  Solution  of  w 
figure  (5)  ; 

title  (  Exact  Solution  for  w') 

plot3  ( coord  ( intma  ,  1 )  ,  coord  ( intma  , 2 )  ,  qex  { 1  }(:),'*' ) 


.1  0.2  0.8;... 
.5  0.3  0.2]; 

16 


92 


1 

3 

5 

7 

9 

11 

13 

15 

17 

19 

21 

23 

25 

27 

29 

31 

33 

35 

37 


o/  o/ 

/o - /o 

%  Function  for  computing  fw  and  fw  with  an  Upwind  Flux  with  p=0 
%  Boundary  Conditions  enforced  on  a  washer  mesh. 

%  Written  by  Benjamin  Davis  and  Asst.  Professor  Jeremy  Kozdon 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 


function  [ fwl , fw2 ,  fw3 ,  fw4 ,  fpl  ,  fp2  ,  fp3  ,  fp4  ]  =  UpwindFlux2DFwFpBC  (Ne,  ql ,  LI 
,  L2 ,  L3 ,  L4 ,  EtoEBC ,  Dx , .  .  . 

Dy,  nx_l  ,  ny_l ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4 ) 

c  =  1; 

for  k  =  l:Ne 

nEl  =  EtoEBC (k,  1 )  ; 
nE2  =  EtoEBC (k, 2)  ; 
nE3  =  EtoEBC(k,3); 
nE4  =  EtoEBC (k, 4)  ; 

%Side  one  of  element  face 
wml  =  Ll*ql  { 1 }  ( k  , : ) 
wpl  =  L3*ql  { 1 }  (nEl , : ) 
if  k  ==  nEl 
wpl  =  -wml ; 

end 

%Side  two  of  element  face 
wm2  =  L2* ql  { 1 } ( k , : ) 
wp2  =  L4*  ql  { 1 }  (nE2  , : ) 
if  k  ==  nE2 
wp2  =  -wm2; 

end 

%Side  three  of  element  face 
wm3  =  L3*ql  { 1 }  (k  , : ) 
wp3  =  Ll*ql  {1 }  (nE3  , : ) 
if  k  ==  nE3 
wp3  =  -wm3; 

end 

%Side  four  of  element  face 
wml  =  L4*ql  { 1 }  (k  , : ) 
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wp4  =  L2*  ql  { 1 }  ( nE4  , : ) 
if  k  ==  nE4 
wp4  =  -wm4; 

end 

%Side  one  of  element  face 

pml  =  nx_l  { k }  *  LI  *Dx{  k }  *  ql  { 2  }  ( k  , : )  '  +  ny_l  { k  [  *  LI  *Dy{  k  [  *  ql  {  2  }  ( k  , : ) 
ppl  =  nx_3 { nEl } *L3*Dx{ nEl } * ql { 2 } ( nEl , : )  '  +  ny_3 { nEl } * L3*Dy { nEl } * ql 
{  2  }  ( nEl  , : ) 

if  k  ==  nEl 
ppl  =  -pml; 

end 

%Side  two  of  element  face 

pm2  =  nx_2  { k }  *  L2  *Dx  {k}*ql{2}(k,:)  '  +  ny_2  { k  }  *  L2  *Dy  { k  }  h-  ql  { 2  }  (  k  , : ) 
pp2  =  nx_4  { nE2  }  *L4*Dx{ nE2  }  *  ql  { 2  }  ( nE2  , : )  '  +  ny_4  { nE2  }  * L4*Dy { nE2  } *  ql 
{ 2  }  ( nE2  , : ) 

if  k  ==  nE2 

pp2  =  -pm2; 

end 

%Side  Three  of  element  face 

pm3  =  nx_3  { k }  *  L3  *Dx  { k }  *  ql  { 2 }  ( k  , : )  '  +  ny_3  { k } *  L3  *Dy  { k  )  *  ql  { 2 }  ( k  , : ) 
pp3  =  nx_l {nE3 } *L1 *Dx{nE3 } * ql { 2 } ( nE3 , : )  '  +  ny_l {nE3 } *  LI *Dy{nE3 } * ql 
{ 2  }  ( nE3  , : ) 

if  k  ==  nE3 

pp3  =  -pm3; 

end 

%Side  Four  of  element  face 

pm4  =  nx_4  { k }  *  L4*Dx{  k }  *  ql  { 2  }  ( k  , : )  '  +  ny_4  { k  [  *  L4*Dy{  k  [  *  ql  {  2  }  ( k  , : ) 
pp4  =  nx_2  { nE4 }  *L2*Dx{  nE4 }  *  ql  { 2  }  ( nE4  , : )  '  +  ny_2  { nE4  }  *  L2*Dy  { nE4  }  *  ql 
{  2  }  ( nE4  , : ) 

if  k  ==  nE4; 
pp4  =  -pm4; 

end 

%Calculation  of  fpl  ,  fp2  ,  fp3  ,  fp4  ,  fwl ,  fw2 ,  fw3 ,  fw4 
f  p  1  { k }  =  (1/2)  *  (pml  -  ppl)  +  l/(2*c)  *(wpl-wml)  ; 

fp2  { k }  =  (1/2)  *(pm2  -  pp2)  +  l/(2*c)  *(wp2-wm2)  ; 

fp3{k[  =  (l/2)*(pm3  -  pp3)  +  l/(2*c)  »(wp3-wm3)  ; 

fp4{k[  =  (l/2)*(pm4  -  pp4)  +  l/(2*c)  *(wp4-wm4)  ; 

fwl{k}  =  (l/2)*(wpl  -  wml)  -  c/2=i(ppl+pml)  ; 
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fw2  { k } 
fw3 1  k } 
fw4  { k  [ 

end 

end 


(l/2)*(wp2  -  wm2) 
(1/2)  *(wp3  -  wm3) 
(l/2)*(wp4  -  wm4) 


-  c/2*(pp2+pm2)  ; 

-  c/2*(pp3+pm3)  ; 

-  c/2*( pp4+pm4)  ; 


O/ 
/O  “ 


o/ 

■  /o 


%Function  for  computing  fw  and  fp  with  and  Upwind  flux  with  no  boundary 
%conditions  being  enforced  .  Used  with  periodic  boundary  conditions  . 
%Written  by  Benjamin  Davis  and  Asst.  Professor  Jeremy  Kozdon 
%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  % 

function  [ fwl , fw2 ,  fw3 ,  fw4 ,  fpl  ,  fp2  ,  fp3  ,  fp4  ]  =  UpwindFlux2DFwFp ( Ne ,  ql  ,L1 , 
L2 ,  L3 ,  L4 ,  EtoE  ,  Dx , .  .  . 

Dy,  nx_l  ,  ny_l ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4 ) 


11 

c  =  1; 

for  k  =  1 

Ne 

13 
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EtoE (k 

,i); 

nE2  = 

EtoE (k 

,2); 

15 

nE3  = 

EtoE (k 

,3)  ; 

nE4  = 

EtoE (k 

,4); 

17 

%Side 
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face 

wml  = 

Ll*ql  { 

1 } ( k , : ) 
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wpl  = 

L3  *  ql { 

1 }(nEl  , : ) 
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/ 
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pml  =  nx_l { k } *  LI *Dx{ k } * ql { 2 } ( k , : )  '  +  ny_l { k } *  LI *Dy{ k } * ql { 2 } ( k , : 
ppl  =  nx_3{nEl }*L3*Dx{nEl }*ql {2}(nEl , : ) '+ny_3{nEl }*L3*Dy{nEl 
[2}(nEl ,:) 

%Side  two  of  element  face 

pm2  =  nx_2 | k } * L2 *Dx{ k } * ql { 2 } ( k  , : )  '  +  ny_2  { k  }  *L2*Dy{  k  }  *  ql  { 2 }  ( k  , : 
pp2  =  nx_4 { nE2 } *L4*Dx{ nE2 } * ql { 2 } ( nE2 , : )  '  +  ny_4 { nE2 } * L4*Dy { nE2 
i 2  }  ( nE2  , : ) 

%Side  Three  of  element  face 

pm3  =  nx_3  { k }  *L3*Dx{  k }  *  ql  { 2  }  ( k  , : )  '  +  ny_3  { k }  * L3*Dy{  k }  *  ql  {  2  }  ( k  , : 
pp3  =  nx_l  { nE3  }  *L1  *Dx{  nE3  }  *  ql  { 2  }  ( nE3  , : )  '  +  ny_l  { nE3  }  *  LI  *Dy  { nE3 
[ 2 }  ( nE3  , : ) 

%Side  Four  of  element  face 

pm4  =  nx_4  { k }  *  L4*Dx{  k }  *  ql  { 2  }  ( k  , : )  '  +  ny_4  { k }  *  L4*Dy{  k }  *  ql  {  2  }  ( k  , : 
pp4  =  nx_2 { nE4 } *L2*Dx{ nE4 } * ql { 2 } ( nE4 , : )  '  +  ny_2 { nE4 } * L2*Dy { nE4 
! 2  }  ( nE4  , : ) 

%Calculation  of  fpl  ,  fp2  ,  fp3  ,  fp4  ,  fwl ,  fw2 ,  fw3 ,  fw4 
f  p  1  { k }  =  ( 1  /  2 )  *  (pml  -  ppl)  +  1/(2*  c)  *(wpl-wml)  ; 

fp2{k[  =  (l/2)*(pm2  -  pp2)  +  l/(2*c)  *(wp2-wm2)  ; 

fp3{k)  =  (l/2)*(pm3  -  pp3)  +  l/(2*c)  *(wp3-wm3)  ; 

fp4  { k }  =  (l/2)*(pm4  -  pp4)  +  l/(2*c)  *(wp4-wm4)  ; 


) 

*  ql 


) 

*  ql 


) 

*  ql 


) 

*  ql 


fwljkj  =  (l/2)*(wpl  -  wml) 

fw2{kj  =  (l/2)*(wp2  -  wm2) 

fw3|k}  =  (l/2)>t(wp3  -  wm3) 

fw4{k[  =  (l/2)>t(wp4  -  wml) 


-  c /2*(ppl+pml)  ; 

-  c/2*(pp2+pm2) ; 

-  c/2*(pp3+pm3)  ; 

-  c/2*( pp4+pm4)  ; 


end 


end 


2 


4 


6 


8 


%Central  Flux  Routine .  Provides  the  center  flux  values  for  fw  and  fp 
%for  the  two-dimensional  problem. 

%Written  by  Benjamin  Davis 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - %  % 
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function  [fwl  ,fw2 ,  fw3 ,  fw4 ,  fpl  ,  fp2  ,  fp3  ,  fp4  ]  =  CenterFlux2DFwFp  (Ne,  ql  ,L1 , 
L2 ,  L3 ,  L4 ,  EtoE  ,  Dx , .  .  . 

Dy,  nx_l  ,  ny_l ,  nx_2  ,  ny_2  ,  nx_3  ,  ny_3  ,  nx_4  ,  ny_4 ) 
for  k  =  l:Ne 

nEl  =  EtoE (k , 1 ) ; 
nE2  =  EtoE (k , 2 ) ; 
nE3  =  EtoE (k ,3) ; 
nE4  =  EtoE(k,4) ; 

%Side  one  of  element  face 
wml  =  LI  *ql  { 1 }  ( k  , : ) 
wpl  =  L3*ql  { 1 }  (nEl  , : ) 

%Side  two  of  element  face 
wm2  =  L2*  ql  { 1 }  ( k  , : ) 
wp2  =  L4*ql  { 1 }  (nE2  , : ) 

%Side  three  of  element  face 
wm3  =  L3*ql  { 1 }  (k  , : ) 
wp3  =  Ll*ql  {1 }  (nE3  , : ) 

%Side  four  of  element  face 
wm4  =  L4*ql  { 1 }  (k  , : ) 
wp4  =  L2*ql  { 1 }  (nE4  , : ) 

fwljk}  =  (l/2)*(wpl  -  wml); 
fw2{k}  =  (1/2)  *(wp2  -  wm2)  ; 
fw3jk|  =  (1/2)  *(wp3  -  wrh3)  ; 
fw4{k[  =  (1/2)  *(wp4  -  wm4)  ; 


%Calculation  of  fpl  ,  fp2  ,  fp3  ,  fp4 
%Side  one  of  element  face 

pml  =  nx_l  {  k  )  x-  LI  *Dx{  k  )  *  ql  { 2  }  ( k  , : )  '  +  ny_l  { k  }  *  LI  *Dy{  k  }  *  ql  {  2  }  ( k  , : ) 
ppl  =  nx_3 { nEl } *L3*Dx{ nEl } * ql { 2 } ( nEl , : )  '  +  ny_3 { nEl } * L3*Dy { nEl } * ql 
! 2 }  ( nEl  , : ) 

%Side  two  of  element  face 

pm2  =  nx_2  { k }  *  L2  *Dx  jk}*ql{2}(k,:)  '  +  ny_2  { k  }  *  L2  *Dy  { k  }  *  ql  {  2  }  ( k  , : ) 
pp2  =  nx_4 { nE2 } *L4*Dx{ nE2 } * ql { 2 } ( nE2 , : )  '  +  ny_4 { nE2 } * L4*Dy{ nE2 } * ql 
I  2  }  ( nE2  , : ) 

%Side  Three  of  element  face 

pm3  =  nx_3  {  k )  *  L3*Dx|  k )  *  ql  { 2  }  ( k  , :  )  '  +  ny_3  { k  }  *  L3*Dy{  k  }  *  ql  {  2  }  ( k  , : ) 
pp3  =  nx_l { nE3 } *L1 »Dx{ nE3 } * ql { 2 } ( nE3 , : )  '  +  ny_l { nE3 } *  LI *Dy  { nE3 } * ql 
1 2  }  ( nE3  , : ) 
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%Side  Four  of  element  face 

pm4  =  nx_4  { k }  *  L4  *Dx  { k }  ql  { 2 }  ( k  , : )  '  +  ny_4  { k }  *  L4  *Dy  { k}*ql{2}(k,:) 
pp4  =  nx_2  { nE4 }  *L2*Dx{  nE4 }  *  ql  { 2 }  ( nE4  , : )  '  +  ny_2  { nE4  }  *  L2*Dy{  nE4  }  *  ql 
i  2  }  ( nE4  , : ) 

fpl  { k }  =  ( 1  / 2 )  *  (pml  -  ppl )  ; 
fp2 { k [  =  (1/2) *(pm2  -  pp2) ; 
fp3 { k [  =  (1/2) *(pm3  -  pp3) ; 
fp4  { k }  =  ( 1  / 2 )  *  (pm4  -  pp4)  ; 

end 

end 
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o/  o/ 

/o - /o 

%  Function  to  solve  the  combined  discretized  equations  from  Chapter  4 
%  for  the  2D  Acoustic  Wave  equation .  Computes  RHS  and  inverts  the 
%  required  matrices  to  compute  w  and  p  for  second  order  wave  equation. 
%  Written  by  Benjamin  Davis  and  Asst.  Professor  Jeremy  Kozdon 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 


function  [R]  =  RHSDG2Dn(Ne,MAT,  Je  ,M,M_1D,M1,  qs  ,Dx,Dy,  LI  ,L2 ,  L3 ,  L4 ,  nx_l , 
ny_l ,  nx_2  ,ny_2  ,  nx_3  ,ny_3  ,nx_4  ,  ny_4  ,  Sjl  ,  S j  2  ,  Sj3  ,  Sj4  ,  fwl ,  fw2 ,  fw3 ,  fw4 , 
fpl  ,  fp2  ,  fp3  ,  fp4  ,  c  , MAT_ksi , MAT_eta ) 

R { 1 }  =  zeros ( size ( qs { 1 })) ; 

R { 2 }  =  zeros ( size ( qs { 2 })) ; 

for  e  =  l:Ne 

R{  2 }  ( e  , : )  =  (MAT{  e }  +M1 '  *  J  e  { e  )  *M)  *  qs  { 1 }  ( e  , : ) 


R{ 2 } ( e , : )  =  R| 
*M_1D*  Sjl  { e  }  *  fwl 

2  }  ( e  , : )  ' 
{e  })  ; 

+  ( ( Dx  { e } 

'* LI ' * nx_l { e}+  Dy{e} 

' *L1 ' * ny_l { e } ) 

R{ 2 } ( e , : )  =  R j 

2  }  ( e  , : )  ' 

+  ( ( Dx  { e } 

' * L2 '* nx_2 { e }+  Dy{e} 

'  *  L2  '  *  ny_2  { e  } ) 

*M_1D*  S j  2  { e  }  h-  fw2  { e  } )  ; 
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R{  2 }  ( e  , : )  =  R{2}(e,:)  '  +  ((Dxje}  L3  '  *  nx_3  {e}+  Dy{e}  '  *L3  '  *  ny_3  { e  } ) 
*M_lD*Sj3  {e}*fw3{e})  ; 

R{  2 }  ( e  , : )  =R{2}(e,:)'  +  ( (Dx{  e  }  ’*  L4  '  *  nx_4  { e}+  Dy{  e  }  '  *L4  ’  *  ny_4  { e  } ) 
*M_1D*  S j  4  { e  }  *  fw4  { e  } )  ; 

R{  1 }  ( e  , : )  =  c  A2*(L1 '  *M_1D*  Sj  1  { e }  ^  f p  1  { e}+  L2  ’  *M_1D*  S j 2  { e  }*  fp2  { e}  +  L3 
'*M_lD*Sj3  {e}*fp3{e}  +  L4'*M_lD*Sj4{e}*fp4{e})  ; 

R{  1 }  ( e  , : )  =  R{  1 )  ( e  , : )  '  -  c  A2*(MAT_ksi{ e}  +  MAT_eta{  e  } )  , qs  { 2 }  ( e  , : ) 

Mil  { e }  =  (MAT)  e }  +M1 '  *  J  e  { e  }  *M)  ; 

[L_1,U_1]  =  lu(Mll{e},  'vector'); 

R1  =  U_1  \  ( L_1  \  R{  2 }  ( e  , : )  '); 

R1  =  R1  ' ; 

R{ 2 } ( e  , : )  =  Rl; 

Ml_2  { e }  =  J  e  { e  }  *M; 

[L_2,U_2]  =  lu(Ml_2{e},  'vector'); 

R2  =  U_2  \  ( L_2  \  R{  1 )  ( e  , : )  ’); 

R2  =  R2  ' ; 

R{  1 }  ( e  , : )  =  R2 ; 

end 

end 


%This  function  computes  the  LGL  grid  and  elements  in  2D.  Modified  by 
3  %Jeremy  Kozdon  and  Benjamin  Davis  for  thesis  project  to  produce  a  skewed 
%mesh . 

s  % 

%Function  given  to  F.X.  Giraldo  '  s  MA4245  class  August  2014. 

7  %Modified  :  Jan  2015 
%Written  by  F.X.  Giraldo  on  4/2008 
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%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  INPUT  LIST :  nelx  and  nely  are  the  number  of  elements  in  x  and  y 
%  nop  is  the  polynomial  order 

%  xgl  are  the  interpolation  points  on  the  element. 

%  OUTPUT  LIST : 

%  coord  are  the  coordinates:  x=coord(:,l)  and  y=coord(:,2) 

%  intma  is  the  connectivity  list  that  points  to  the  global 

%  gridpoint  number 

%  bsido  is  the  boundary  data  (used  by  ISIDE  and  FACE) 

%  iperiodic  points  to  another  point  if  periodicity  is 

%  applicable 

%  npoin  =  number  of  global  points 

%  nelem  =  number  of  elements 

%  nboun  =  number  of  boundary  edges 

%  nface  =  number  of  faces /edges  in  the  grid 

Of  0/ 

/o - /o 


function  [coord  ,  intma  ,  bsido  ,  iperiodic  ,  npoin  ,  nelem  ,  nboun ,  nface  ]  = 
create_grid_2d(  nelx  ,  nely  ,nop ,  xgl  ,  plot_grid  ,  skew_mesh) 
%Define  Grid  Dimensions 
ngl=nop  +  l; 

npoin  =  (nop*nelx  +  1)  *(  nop  *  nely  +  1); 
nelem=nelx  *  nely ; 
nboun=2*nelx  +  2*nely; 
nface  =2*nelem  +  nelx  +  nely; 

%Initialize  Global  Arrays 
coord=zeros (npoin  ,2) ; 
intma=zeros  (nelem  ,  ngl  ,  ngl )  ; 
bsido=zeros (nboun, 4) ; 
iperiodic=zeros  (npoin  ,  1 )  ; 

%Initialize  Local  Arrays 
node=zeros (npoin  ,  npoin)  ; 

%Set  some  constants 
xmin  =  - 1 ; 
xmax=+l; 
ymin  =  - 1 ; 
ymax=+l; 

dx  =  (xmax-xmin)  /  nelx  ; 
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dy=(ymax-ymin)  /  nely ; 
nop=ngl - 1 ; 
nx=nelx  *nop  +  1; 
ny=nely  *nop  +  1; 

°/(GENERATE  COCKD 
ip  -0; 

jj  =0; 

for  k  =  l:nely 

y0=ymin  +  real (k-1) *dy; 
if  (k  ==  1) 

11  =1; 
else 

11  =2; 

end 

for  1=11 : ngl 

i)=jj  +i; 

i  i  =0; 

for  i=l:nelx 

x0=xmin  +  real(i-l)*dx; 
xc  =  xO  +  [0,l;0,l]*dx; 
yc  =  yO  +  [0  ,0;  1 
if  (skew_mesh=  =  l) 

xc=xc  +  ( 1/8)  *  sin  (  pi  *yc  )  .  *(  1  -  xc  )  .*(  1  +  xc  )  ;%x ; 
yc=yc  +( 1/8)  *  sin  ( pi  *xc )  .*(  1  -  yc )  .*(  1  +  yc )  ;%y ; 
end 

if  (i  ==  1) 
jl=l; 
else 
jl=2; 

end 

for  j  =j 1  : ngl 
i  i  =  i  i  +  1 ; 
ip=ip  +  1; 
ax=(xgl(j  )+l)/2; 

ay =( xgl ( 1 ) +l)/2; 

x  =  xc(l,l)*(l-ax)*(l-ay)  +  xc(l,2)*ax*(l 
+  xc(2,l)*(l-ax)*(  ay)  +  xc(2,2)*ax*( 
y  =  yc(l,l)*(l-ax)*(l-ay)  +  yc(l,2)*ax=f(l 


ay)  . . . 

ay) ; 

ay)  .  . . 
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+  yc(2 ,1) *(1 -ax) *(  ay)  +  yc (2 ,2) *ax*(  ay); 
coord (ip  , 1 ) =x; 
coord (ip  ,2) =y ; 
if  (skew_mesh==2) 

coord(ip,l)=x  +  (l/8)*sin(pi*y)*(l-x)*(l  +  x)  ;%x ; 
coord(ip,2)=y  +  (l/8)*sin(pi*x)*(l-y)*(l+y)  ;%y ; 
end 

node( ii  , j j )  =  ip ; 
end  %j 

end  %i 
end  %1 
end  %k 

“/(GENERATE  MMA 
ie  =0; 

for  k  =  l:nely 
for  i =1 : nelx 
ie  =  ie  +1; 
for  l=l:ngl 

jj  =(ngl -1) *(k-l)  +  1 ; 

for  j  =1 :  ngl 

ii=(ngl-l)  *(i  -1)  +  j  ; 
ip=node(  ii  ,  j  j  ) ; 
intma(  ie  ,  j  ,1  )=ip  ; 
end  %j 
end  %1 
end  %i 
end  %k 

%Generate  BSIDO 
ib  =0; 

for  i=l:nelx 
ie  =  i ; 
ib  =  ib  +1; 

il=(i-l)*(  ngl  - 1 )  +  1; 
i2=(i  -1)  *(ngl  -1)  +  ngl; 
ipl=node( il  ,1 ) ; 
ip2=node(  i2  ,1 ) ; 
bsido  ( ib  ,  1 )  =ipl  ; 
bsido ( ib  ,2)  =ip2  ; 
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bsido  ( ib  ,3 )  =ie  ; 
bsido  ( ib  ,4)  =  6; 

end 

%Right  Boundary 
for  i=l:nely 

ie  =(nelx  )  *(  i ) ; 
ib=ib  +1; 

il  =(i  -1)  *(ngl  -1)  +  1; 
i2=(i  -1)  *  ( ngl  - 1 )  +  ngl ; 
ipl=node(nx, il  ) ; 
ip2=node(nx,  i2  ) ; 
bsido  ( ib  ,  1 )  =ipl  ; 
bsido ( ib  ,2)  =ip2  ; 
bsido ( ib  ,3 )  =ie  ; 
bsido  ( ib  ,4)  =6; 

end 

%Top  Boundary 
for  i=nelx  :  -1 : 1 

ie=nelem  -  (nelx  -  i); 
ib  =  ib  +1; 

11  =(i  -1)  *(ngl  -1)  +  ngl; 

1 2  = ( i  -1) * ( ngl - 1 )  +  1; 
ipl=node(il  , ny) ; 
ip2=node(i2  ,ny) ; 
bsido  ( ib  ,  1 )  =ipl  ; 
bsido  ( ib  ,2)  =ip2  ; 
bsido  ( ib  ,3 )  =ie  ; 
bsido  ( ib  ,4)  =6; 

end 

%Left  Boundary 
for  i=nely : -1 : 1 

ie =(nelx ) *( i -1 )  +  1; 
ib  =  ib  +1; 

11  =(i  -1)  *(ngl  -1)  +  ngl; 

1 2  = ( i  -1) *  ( ngl - 1 )  +  1; 
ipl=node(l , il ) ; 
ip2=node(l , i2 ) ; 
bsido  ( ib  ,  1 )  =ipl  ; 
bsido ( ib  ,2)  =ip2  ; 
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bsido  ( ib  ,3 )  =ie  ; 
bsido  ( ib  ,4)  =6; 

end 

%Perio  dicity 
for  i=l:npoin 

iperiodic  ( i  )  =  i ; 

end 

%X-  Periodicity 
for  i=l:ny 

il=node (1 , i ) ; 
i2=node(nx, i ) ; 
iperiodic  ( i 2  )  =  i  1  ; 

end 

%Y-  Periodicity 
for  i=l:nx 

il =node ( i , 1 ) ; 
i2=node ( i ,ny ) ; 

iperiodic(i2)  =  iperiodic(il  )  ; 

end 

%Plot  Grid 
if  (plot_grid  ==  1) 
x=zeros (5,1) ; 
y=zeros (5,1)  ; 
figure  ; 
hold  on; 
for  e  =  l:nelem 

for  j  =l:ngl  -1 

for  i=l:ngl-l 

il  =intma  ( e  ,  i  ,  j  )  ; 
i2=intma  ( e  ,  i  +1 ,  j  )  ; 
i3=intma (e,  i  +1, j  +1) ; 
i4=intma  ( e  ,  i  ,  j  +1) ; 

x( 1 ) =coord ( il ,1) ;  y ( 1 ) =coord ( il ,2) ; 
x(2) =coord ( i2 ,1) ;  y(2)=coord ( i2 ,2) ; 
x(3) =coord ( i3 ,1) ;  y(3)=coord ( i3 ,2) ; 
x(4) =coord ( i4 , 1 )  ;  y (4 ) =coord ( i4  ,2 )  ; 
x(5) =coord ( il ,1) ;  y(5)=coord ( il ,2) ; 
plot_handle=plot (x,y, ' -r ' ) ; 
set  ( plot_handle  ,  '  LineWidth  '  ,1.5)  ; 
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end 

end 

i  1  =intma  ( e  ,  1  ,  1 )  ; 
i2=intma  (e  ,  ngl  ,  1 )  ; 
i3=intma  (e  ,  ngl  ,  ngl )  ; 
i4=intma(e,l  ,ngl)  ; 

x( 1 ) =coord ( il  ,1) ;  y ( 1 ) =coord ( il  ,2)  ; 
x(2) =coord  ( i2  ,1)  ;  y (2 ) =coord ( i2  ,2)  ; 
x(3) =coord ( i3  ,1)  ;  y (3) =coord ( i3 ,2)  ; 
x(4)=coord(i4 ,1)  ;  y (4) =coord ( i4 ,2 ) ; 
x(5) =coord ( il  ,1) ;  y (5 ) =coord ( il  ,2)  ; 
plot_handle  =  plot  (x,y,  '  -b  '  )  ; 
set  ( plot_handle  ,  '  LineWidth  ,2)  ; 

end 

title_text  =[  '  Grid  Plot  For:  Ne  =  '  num2str  (nelem)  N=  '  num2str( 
nop )  ] ; 

title  ([  title_text  ]  ,  'FontSize  '  ,18) ; 
xlabel(  'X'  ,  'FontSize  '  ,18)  ; 
ylabel(  'Y'  ,  'FontSize  '  ,18)  ; 
axis  image 

end 
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%Function  provided  to  Professor  F.X.  Giraldo  '  s  MA4245  class. 

%Used  by  Benjamin  Davis  for  his  thesis  project  . 

%This  subroutine  creates  the  array  ISIDE  which  stores  all  of 
%the  information  concerning  the  sides  of  all  the  elements. 

%Written  by  Francis  X.  Giraldo  on  1/01 
%  Naval  Postgraduate  School 

%  Department  of  Applied  Mathematics 

%  Monterey,  CA  93943-5502 

%  INPUT  LIST :  intma  =  element  connectivity 

%  bsido  =  boundary  info  (which  points  are  on  a  boundary, 

%  which  element  it  belongs  to  and  the  boundary 

%  condition )  . 

%  npoin  =  number  of  global  points 
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%  nelem  =  number  of  elements 

%  nboun  =  number  of  boundary  faces /edges 

%  nside=nface  are  the  number  of  sides / face /edges  in  the  grid 

%  ngl  =  number  of  points  in  one  direction  in  an  element 

%  OUTPUT  LIST :  iside  =  face  information  such  as  which  points  are  on  a 
%  face  which  elements  they  belong  to  and,  if  a 

%  boundary,  what  is  the  boundary  condition. 

%  jeside  =  for  each  element  and  each  edge  gives  the  FACE 

%  number 

o/  o/ 

/o - /o 

function  [  iside  ,  j  eside  ]  =  create_side  (intma  ,  bsido  ,  npoin  , nelem  , nboun. 


nface  ,  ngl ) 

%global  arrays 
iside  =  zeros (nface  ,4)  ; 
jeside=  zeros  (nelem  ,4)  ; 

%local  arrays 
lwher  =  zeros  ( npoin  ,  1 )  ; 
lhowm  =  zeros  ( npoin  ,  1 )  ; 
icone  =  zeros  (5*npoin  ,  1 )  ; 
inode  =  zeros (4,1); 
jnode  =  zeros (4,1)  ; 

%Fix  Inode 
inode  ( 1 )  =1; 
inode  (2)=ngl ; 
inode  (3)=ngl ; 
inode  (4)  =  1; 
jnode  ( 1 )  =  1; 
jnode  (2)  =  1; 
jnode  (3)=ngl ; 
jnode  (4)  =ngl ; 

%count  how  many  elements  own  each  node 
for  in=l:4 

for  ie=l:nelem 

ip=intma  ( ie  ,  inode  ( in  )  ,  jnode  ( in  ) )  ; 
lhowm(  ip  )=lhowm(  ip  )  +  1; 
end  %ie 
end  %in 

%track  elements  owning  each  node 
lwher  ( 1 )  =0; 
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for  ip=2:npoin 

lwher  ( ip  )=lwher  ( ip  - 1 )  +  lhowm(ip-l); 
end  %ip 

%another  tracker  array 
lhowm  =  zeros  ( npoin  ,  1 )  ; 
for  in  =  l:4 

for  ie  =1  :nelem 

ip=intma  ( ie  ,  inode  ( in  )  ,  jnode  ( in  ) )  ; 
lhowm  (ip)  =lhowm  ( ip  )  +  1 ; 
jloca=lwher  ( ip  )  +  lhowm  (ip); 
icone ( j  loca )=ie ; 
end  %ie 
end  %in 

°/<LOOP  OVER  THE  NODES 
iloca  =0; 
for  ip=l:npoin 
ilocl  =  iloca  ; 
i  e  1  e  =lhowm ( ip  )  ; 
if  ( iele  ~=  0  ) 

iwher=lwher  ( ip  )  ; 

“/LOOP  OVER  THOSE  ELEMENTS  SURROUNDING  NODE  IP 
ip  1  =  ip  ; 

for  iel  =1:  iele 

ie  =  icone  ( iwher+iel )  ; 

%find  out  position  of  ip  in  intma 
for  in  =  l:4 
inl=in ; 

ipt=intma  ( ie  ,  inode  ( in  )  ,  jnode  ( in  )  )  ; 
if  (ipt  ==  ip) 
break 

end 
end  %in 

%Check  Edge  of  Element  IE  which  claims  IP 
1=0; 

for  jnod  =  l:2:3 
iold  =0; 

j=j +1; 

in2=in  +  jnod ; 
if  (in2  >  4) 
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in2=in2 -4; 


end 

ip2=intma  ( ie  ,  inode  ( in2  )  ,  jnode  ( in2  )  )  ; 
if  ( ip2  >=  ipl ) 

%check  whether  side  is  old  or  new 
if  ( iloca  ~=  ilocl  ) 

for  is  =  ilocl +1:  iloca 
iside  (is  ,2)  ; 
j  1  o  c  a  =  i  s  ; 

if  (iside(is,2)  ==  ip2 ) 
iold  =1; 
break ; 

end 
end  %is 
end  %iloca 
if  (iold  ==  0) 

°XNEW  SIDE 
iloca  =  iloca  +  1; 
iside(iloca  ,l)=ipl  ; 
iside(iloca  ,2)=ip2; 
iside  (iloca  ,2+j  )  =  ie  ; 
elseif  ( iold  ==  1) 

°/<DLD  SIDE 

iside(jloca  ,2+j  )  =  ie  ; 
end  %iold 
end  %ip2 
end  %jnod 
end  %i e  1 

%Perform  some  Shifting  to  order  the  nodes  of  a  side  in  COV 
direction 

for  is  =  ilocl  +1:  iloca 

if  ( iside  ( is  ,3)  ==  0) 

iside(is  ,3)=iside(is  ,4)  ; 
iside  ( is  ,4)  =0; 
iside(is  ,l)=iside(is  ,2)  ; 
iside  ( is  ,2)=ipl ; 
end  %iside 
end  %is 
end  %if  iele 
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end  %ip 

if  (iloca  ~=  nface) 

disp  (  'Error  in  SIDE,  iloca  nface  =  '); 

iloca 

nface 

pause 

end 

%RESET  THE  BOUNDARY  MARKERS 
for  is=l:nface 

if  ( iside  ( is  ,4)  ==  0) 
i  1  =iside  ( is  ,  1 )  ; 
ir  =  iside  (is  ,2)  ; 
ie  =  iside  (is  ,3)  ; 
for  ib=l:nboun 

ibe=bsido  (ib  ,3)  ; 
ibc=bsido  ( ib  ,4)  ; 
if  (ibe  ==  ie) 

ilb=bsido  ( ib  ,  1 )  ; 
irb=bsido  ( ib  ,2  )  ; 
if  (  ilb  ==  il  &&  irb  ==  ir  ) 
iside  (is  ,4)  =-ibc  ; 
break 
end  %i  1  b 
end  %ibe 
end  %ib 
end  %iside 
end  %is 

°/<K)RM  ELEMENT/ SIDE  CONNECTIVITY  ARRAY 
for  is  =1:  nface 

iel=iside  (is  ,3)  ; 
ier  =  iside  (is  ,4)  ; 
is  1  =iside  (is  ,  1 )  ; 
is2  =  iside  (is  ,2)  ; 

%LEFT  SIDE 
for  in=l:4 

il  =intma  ( iel  ,  inode  ( in  )  ,  jnode  ( in  ) )  ; 
inl=in  +  1; 
if  (ini  >  4) 
ini  =1; 
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end  %inl 

i2=intma  ( iel  ,  inode  ( ini )  ,  jnode  ( ini ) )  ; 
if  ( ( i  s  1  ==  il)  &&  (is2  ==  i2  ) ) 
jeside  ( iel  ,in)  =  is; 
end  %i  s 1 
end  %in 
%RIGHT  SIDE 
if  (ier  >  0) 
for  in=l:4 

i  1  =intma  ( ier  ,  inode  ( in  )  ,  jnode  ( in  )  )  ; 
inl=in  +  1; 
if  (ini  >  4) 
ini  =1; 
end  %inl 

i2=intma  ( ier  ,  inode  ( ini  )  ,  jnode  ( ini  )  )  ; 
if  ( (  is  1  ==  i2 )  &&  ( is2  ==  il)) 
jeside  ( ier  ,in)  =  is  ; 
end  %i  s 1 
end  %in 
end  %ier 
end  %is 


i 


3 


5 


7 


9 


11 


13 


%Function  provided  to  Professor  F.X.  Giraldo  's  MA4245  class.  Used  by 
%Benjamin  Davis .  This  subroutine  constructs  the  Side  Information  for 
%a  High  Order  Spectal  Element  Quads 
%Written  by  Francis  X.  Giraldo 
%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  INPUT  LIST :  iside  =  face  information  to  know  which  points  are  on  a 
%  face  and  which  elements  they  belong  to 

%  intma  =  element  connectivity 

%  nface  =  number  of  faces /edges  in  the  grid 

%  ngl  =  number  of  interpolation  points  in  one  direction 

%  in  an  element 
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%  OUTPUT  LIST :  face  =  face  information  stating  which  local  edge  number 


%  the  face  is  on  and  which  elements  they  belong  to 

%  (left  and  right  neighbors ).  are  the  metric  terms 

%  needed  to  imapl  and  imapr  give  the  tensor  -  product  (i,j) 

%  points  of  the  edge  points  . 

0/  0/ 

/o - /o 


function  [face,  imapl ,  imapr  ]=  create_face(iside  ,  intma  ,  nface  ,  ngl ) 

%global  arrays 

face  =  zeros (nface ,4)  ; 

imapl=zeros  (4,2, ngl )  ; 

imapr =zer os  (4  ,2  ,  ngl )  ; 

%local  arrays 
inode=zeros (4,1)  ; 
jnode=zeros (4,1) ; 

%Construct  Boundary  Pointer 

inode  ( 1 )  =1; 

inode  (2)=ngl ; 

inode  (3)=ngl ; 

inode  (4)  =  1; 

jnode  ( 1 )  =1; 

jnode  (2)  =1; 

jnode  (3)=ngl ; 

jnode  (4)  =ngl ; 

%Construct  IMAP  arrays 
for  l=l:ngl 
%eta  =-l 

imapl  (1  ,1  , 1  )  =  1  ; 
imapl  (1  ,2 , 1 )  =1; 
imapr  (1,1,1)  =ngl  +  1  - 1  ; 
imapr  (1  ,2 , 1 )  =1; 

%ksi=+l 

imapl  (2 ,1  , 1  )=ngl ; 
imapl  (2  ,2 , 1  )  =  1  ; 
imapr  (2,1  , 1  )=ngl; 
imapr  (2,2,1  )=ngl  +  l- 1  ; 

%eta=+l 

imapl  (3,1,1  )=ngl  +  l- 1  ; 
imapl  (3,2,1  )=ngl ; 
imapr  (3  ,1,1  )  =  1 ; 
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imapr  (3  ,2 , 1  )=ngl; 

%ksi  =-l 

imapl  (4,1  , 1 )  =1; 
imapl  (4,2,1)  =ngl  +  1  - 1  ; 
imapr  (4  ,1  , 1 )  =1; 
imapr  (4  ,2 , 1  )  =  1  ; 
end  %1 

%loop  thru  the  sides 
for  i=l:nface 

ipl  =  iside  ( i  ,1 )  ; 
ip2  =  iside  ( i  ,2 )  ; 
iel =iside ( i ,3 ) ; 
ier  =  iside  ( i  ,4)  ; 

%check  for  position  on  Left  Element 
for  j  =1:4 
jl=j ; 
j  2 = j  +1; 
if  (j2  >  4) 

i2  =i; 

end  %j2 

jpl=intma  ( iel  ,  inode  ( j  1  )  ,  jnode  ( j  1  ) )  ; 
jp2=intma  ( iel  ,  inode  ( j2  )  , jnode  ( j2  ) )  ; 
if  (ipl  ==  jpl  &&  ip2  ==  jp2) 
face  (i  ,  1 )  =  j  ; 
break ;  %leave  J  loop 
end  %ipl 
end  %j 

%check  for  position  on  Right  Element 
if  (ier  >  0) 
for  j  =1 :4 
jl=j ; 
j2=j+l; 
if  ()2  >  4) 

]2=1; 
end  %j2 

jpl=intma  ( ier  ,  inode  ( j  1  )  ,  jnode  ( j  1  )  )  ; 
jp2=intma  ( ier  ,  inode  ( j2  )  ,  jnode  ( j2  )  )  ; 

if  (ipl  ==  j p 2  &&  ip2  ==  jpl) 
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%Store  Elements 
face ( i ,3)=iel ; 
face(i  ,4)  =  ier ; 
end  %i 


into 


face 


0/ 
/O  “ 


O/ 

■  /o 


%Function  for  changing  the  face  code  to  enforce  p=0  boundary  conditions 
%for  2nd  order  acoustic  wave  equation  on  a  washer. 

%Written  by  Benjamin  Davis  Created  :  09  May  2015 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function  [  face  ]  =  faceBC  ( nface  ,  f ace  ) 
for  k  =  1:  nface 

if  face (k, 1 )  ==  2 

if  face(k,4)  ==  -6 
face(k,4)  =  -1; 

end 

end 

end 

for  k  =  1:  nface 

if  face (k, 1 )  ==  4 

if  face(k,4)  ==  -6 
face  (k,4)  =  -1; 

end 

end 

end 

end 
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0/  0/ 

/o - /o 

%Function  provided  by  Professor  F.X.  Giraldo  to  his  MA4245  class. 

%Used  by  Benjamin  Davis 

%This  subroutine  builds  Periodic  BCs  along  the  4  edges  of  a  rectangular 
%domain . 

%Written  by  Francis  X.  Giraldo  on  2/2007 


%  Naval  Postgraduate  School 

%  Department  of  Applied  Mathematics 

%  Monterey,  CA  93943-5216 

%  INPUT  LIST :  iside  =  face  information 
%  face  =  more  face  information 

%  coord  =  gridpoint  coordinates 

%  nface  =  number  of  faces 

%  nboun  =  number  of  boundaries 

%  OUTPUT  LIST :  face  =  augments  the  FACE  data  structure  to  include 

%  periodicity 

O/  O/ 

/o - /o 


function  face  =  create_face_periodicity  ( iside  ,  face  ,  coord  ,  nface  , nboun) 

%Constant 

tol=le-6; 

%Local  arrays 
ilef t=zeros (nboun/4 ,1 ) ; 
iright=zeros (nboun/4 ,1 ) ; 
itop  =  zeros  (nboun/4 , 1 )  ; 
ibot=zeros (nboun/4 , 1 ) ; 

%initialize 

nleft=0;  nright=0;  ntop=0;  nbot=0; 

%Find  Extrema  of  Domain 

xmax=max(  coord  ( :  ,  1 )  )  ;  xmin=min(  coord  ( :  ,  1 )  )  ; 
ymax=max(  coord  ( :  ,  2  )  )  ;  ymin=min(  coord  ( :  ,  2  )  )  ; 

%loop  thru  sides  and  extract  Left,  Right,  Bot ,  and  Top 
for  is  =1:  nface 

“/Check  for  Periodicity  Edges 
ier  =  face (is  ,4 )  ; 
if  (ier  ==  -6) 

il  =  iside  (is  ,1)  ;  i2  =  iside  ( is  ,2)  ; 
xm=0.5*(  coord(il,l)  +  coord(i2,l)  ); 
ym=0.5*(  coord(il,2)  +  coord(i2,2)  ); 


114 


40 

42 

44 

46 

48 

50 

52 

54 

56 

58 

60 

62 

64 

66 

68 

70 

72 

74 

76 


%check  Grid  Point 

if  (  abs  (xm  -  xmin)  <  tol  )  %left  boundary 
nleft  =  nleft  +  1; 
ilef t (  nleft )  =  is ; 

elseif  (  abs(xm  -  xmax)  <  tol  )  %right  boundary 
nright=nright  +  1; 
iright(nright)  =  is  ; 

elseif  (  abs(ym  -  ymin)  <  tol  )  %bottom  boundary 
nbot=nbot  +  1 ; 
ibot (nbot)=is ; 

elseif  (  abs(ym  -  ymax)  <  tol  )  %top  boundary 
ntop=ntop  +  1; 
itop  (ntop  )  =  is  ; 
else 

disp  (  'No  match  in  PERIODIC_BCS  for  is  ier  =  '); 
is 
ier 
pause 
end  %if 
end  %ier 
end  %is 

%Loop  through  Periodic  BCs 
%First  :  Do  Left  and  Right 
for  i=l:nleft 

isl=ilef t (i) ; 
il  =iside  ( isl  ,  1 )  ; 
yll=coord  ( il  ,2)  ; 

%Search  for  Corresponding  Right  Edge 
for  j  =1 :  nright 

isr  =  iright  ( j  )  ; 
i2  =  iside  ( isr  ,2)  ; 
yr2=coord ( i2 ,2)  ; 

if  (  abs(yll-yr2)  <  tol  )  %they  match 
face(isl  ,2)=face(isr  ,1)  ; 
face(isl  ,4)=face(isr  ,3)  ; 

face  ( isr  ,3)  =  -6;  %means  skip  it  due  to  Periodicity 
iside(isl  ,4)=iside(isr  ,3)  ; 
iside  ( isr  ,3)  =  -6; 
break ; 
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end  %if 
end  %j 
end  %i 

%Second :  Do  Top  and  Bottom 
for  i=l:nbot 

i s  1  =ibot  ( i ) ; 
il  =iside  ( isl  ,  1 )  ; 
xll=coord  ( il  ,  1 )  ; 

%Search  for  Corresponding  Top  Edge 
for  j  =l:ntop 

isr=itop  ( j  )  ; 
i2  =  iside  ( isr  ,2)  ; 
xr2=coord ( i2  , 1 )  ; 

if  (  abs(xll-xr2)  <  tol  )  %they  match 
face(isl  ,2)=face(isr  ,1)  ; 
face(isl  ,4)=face(isr  ,3)  ; 

face  ( isr  ,3)  =  -6;  %means  skip  it  due  to  Periodicity 
iside  ( isl  ,4)=iside  ( isr  ,3)  ; 
iside  (isr  ,3)  =  -6; 
break ; 
end  %if 
end  %j 
end  %i 
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%Function  given  to  Professor  F.X.  Giraldo  '  s  MA4245  class.  Modified  by 
%Ben  Davis  and  Jeremy  Kozdon. 

%This  function  computes  the  Metric  Terms  for  General  2D  Quad  Grids . 
%Written  by  F.X.  Giraldo  on  4/2008 
%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  INPUT  LIST :  coord  =  gridpoint  coordinates 
%  intma  =  element  connectivity 

%  psi  =  basis  functions  defined  as  psi(NGL,NQ) 

%  dpsi  =  derivative  of  basis  functions  defined  as 
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%  dpsi  (NGL,NQ) 

%  nelem  =  number  of  elements 

%  ngl  =  number  of  interpolation  points  in  one  direction  in 

%  an  element 

%  nq  =  number  of  integration / quadrature  points  in  one 

%  direction  in  an  element. 

%  OUTPUT  LIST :  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  =  are  the  metric  terms  needed  to 
%  compute  derivatives  in 

%  physical  space 

%  x_ksi  ,  x_eta  ,  y_ksi  ,  y_eta  =  are  the  metric  terms  needed  to 

%  compute  derivatives  in 

%  physical  space 

%  jac  =  weights  *  Jacobian  defined  as  j  ac  (nelem  ,nq,nq) 

O/  0/ 

/o - /o 

function  [  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  ,  x_ksi  ,  x_eta  ,  y_ksi  ,  y_eta  ,  jac  ]  = 


metrics2(  coord  ,  intma  ,  psi  ,  dpsi  ,  nelem  ,  ngl ,  nq) 
%Initialize  Global  Arrays 
ksi_x=zeros  (nelem  ,nq,nq)  ; 
ksi_y=zeros  (nelem  ,nq,nq)  ; 
eta_x=zeros  (nelem  ,nq,nq)  ; 
eta_y=zeros  (nelem  ,nq,nq)  ; 
jac=zeros  (nelem  ,nq,nq)  ; 

%Initialize  Local  Arrays 
x_ksi=zeros  (nelem  ,nq,nq)  ; 
x_eta=zeros  (nelem  ,nq,nq)  ; 
y_ksi=zeros  (nelem  ,nq,nq)  ; 
y_eta=zeros  (nelem  ,nq,nq)  ; 
x=zeros ( ngl , ngl )  ; 
y=zeros ( ngl ,ngl)  ; 

%loop  thru  the  elements 
for  ie=l:nelem 

%Store  Element  Variables 
for  j  =l:ngl 
for  i=l:ngl 

ip=intma  ( ie  ,  i  ,  j  )  ; 
x(i , j ) =coord ( ip  ,1)  ; 
y(i , j ) =coord ( ip  ,2)  ; 
end  %i 
end  %j 
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%Construct  Mapping  Derivatives:  dx/dksi  ,  dx/deta  ,  dy/dksi,  dy/deta 
[x_ksi(ie  ,:)  ,x_eta(ie  ,:)  ]  =  map_deriv  ( psi  ,  dpsi  ,x,ngl  ,  nq)  ; 

[y_ksi  ( ie  , :  ,:)  ,y_eta  ( ie  , :  , : )  ]  =  map_deriv  ( psi  ,  dpsi  ,y,ngl  ,nq)  ; 
%Construct  Inverse  Mapping:  dksi/dx,  dksi/dy,  deta/dx,  deta/dy 
for  j  =l:nq 
for  i =1  :nq 

xjac  =  x_ksi  ( ie  ,  i  ,  j  )  *y_eta  ( ie  ,  i  ,  j  )  -  x_eta  ( ie  ,  i  ,  j  )  *y_ksi  ( ie  ,  i  ,  j  ) ; 
ksi_x(ie  ,  i  ,  j  )=  +  1.0/xjac*y_eta(ie  ,  i  ,  j  )  ; 
ksi_y(ie  ,  i  ,  j  )  =  -1.0/xjac*x_eta(ie  ,  i  ,  j  )  ; 
eta_x(ie  ,i  ,  j  )  =  -1.0/xjac  *y_ksi  (ie  ,  i  ,  j  ) ; 
eta_y(ie  ,  i  ,j  )=  +  1.0/xjac*x_ksi(ie  ,  i  ,  j  )  ; 
jac (ie ,i , j )=abs(xjac) ; 
end  %i 
end  %j 
end  %ie 
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%  Function  created  for  taking  the  x  coordinates  from  coord  and 
%  turning  them  into  radius  polar  coordinates  based  on  user  inputs  . 
%Written  by  Benjamin  Davis  %Created :  14  April  2015 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function[R]  =  radius  ( coord  ,  rl  ,  r2  ) 
n  =  length ( coord (:, 1 )) ; 
for  i  =  l:n 

R(i,l)  =  ( ( r2  -  rl  ) /2 ) *coord ( i  ,  1 )  +  (rl  +  r 2 ) / 2 ; 

end 

end 
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/o - /o 

%Function  to  turn  the  y  coord  into  theta  .  Conversion  from  Cartesian  to 
%polar  coordinates.  Also,  rn  is  a  scaling  portion  of  the  washer. 
%Written  by  Benjamin  Davis  Created :  14  April  2015 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function  [theta]  =  thetapolar  ( coord  ,  rn ) 
n  =  length ( coord  (:  ,2 ))  ; 

rl  =  0; 

r2  =  (2*pi)/rn; 
for  i  =  l:n 


theta(i,l)  =  ( ( r2 - r 1 ) /2 ) *  coord ( i , 2 )  +  (rl  +  r 2 ) / 2 ; 

end 

end 
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O/  0/ 

/o - /o 

%Function  to  overwrite  original  coord  for  the  conversion  to  polar  mesh. 
%Function  is  required  for  the  washer  mesh. 

%Written  by  Benjamin  Davis  Created :  14  April  2015 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 

function  [  coord  ]  =  newcoord  (R,  theta  ,  coord  ) 
n  =  length ( coord (:, 1 )) ; 
for  i  =  l:n 

coord(i,l)  =  R( i ) *cos ( theta ( i )) ; 
coord(i,2)  =  R(  i )  *  sin ( theta ( i )) ; 

end 

end 
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o/  o/ 

/o - /o 

%Element  to  Element  function  with  Boundary  Conditions  for  Washer  Grid 
%Synopsis  :  Element  to  element  function  for  identifying  and  storing 
%information  for  looking  across  the  faces  of  an  element  to  it's 
%corresponding  neighbor  . 

%Written  by  Benjamin  Davis 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%  pL  =  Face  of  Left  Element 

%  pR  =  Face  of  Right  Element 

%  Ls  =  Number  of  Left  Element 

%  Rs  =  Number  of  Right  Element 

%Output :  EtoE  is  a  Matrix  of  Number  of  Elements  (Row)  vs.  the  Number 

%of  Sides  per  element  (Four  for  a  square  grid).  For  each  element,  it 
%displays  its  corresponding  neighbor  with  respect  to  boundary 
%conditions  . 

O/  O/ 

/o - /o 


function  [EtoEBC]  =  EtoEfunctBC ( nface  ,  f ace  ,Ne) 
EtoEBC  =  zeros  (Ne, 4); 
for  1=1:  nface 

pL  =  face(l,l);  %Face  of  left  element 
pR  =  face(l,2)  ;  %Face  of  Right  Element 
Ls  =  face(l,3)  ;  %Number  of  Left  Element 
Rs  =  face(l,4);  %Number  of  Right  Element 

if  RS  ==  -1 

EtoEBC  ( Ls  ,pL)  =  Ls  ; 
elseif  (Ls  ~=  -6  &&  Rs  ~=  -6) 

EtoEBC  ( Ls  ,pL)  =  Rs ; 

EtoEBC  (Rs  ,pR)  =  Ls  ; 

end 

end 

end 
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o/  o/ 

/o - /o 

%Element  to  Element  function  for  a  square  grid  with  periodic  boundary 
%conditions  . 

°/o/oSynopsis  :  Element  to  element  function  for  identifying  and  storing 
%information  for  looking  across  the  faces  of  an  element  to  it  '  s 
%corresponding  neighbor  . 

%Written  by  Benjamin  Davis 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%  pL  =  Face  of  Left  Element 

%  pR  =  Face  of  Right  Element 

%  Ls  =  Number  of  Left  Element 

%  Rs  =  Number  of  Right  Element 

O/ 

/O 

%Output :  EtoE  is  a  Matrix  of  of  Number  of  Elements  (Row)  vs.  the 

Number 

%of  Sides  per  element  (Four  for  a  square  grid).  For  each  element,  it 

%displays  its  corresponding  neighbor. 

0/  0/ 

/o - /o 

function  [EtoE]  =  EtoEfunct  ( nface  ,  face  ,Ne) 

EtoE  =  zeros  (Ne, 4); 
for  1=1:  nface 


pL  =  face (1,1)  ; 
pR  =  face (1,2)  ; 

Ls  =  face (1,3)  ; 

Rs  =  face (1 ,4)  ; 

if  ( Ls  ~=  -6  &&  Rs  ~=  -6) 
EtoE(Ls,pL)  =  Rs ; 
EtoE(Rs,pR)  =  Ls ; 

end 

end 

end 
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%Function  for  building  the  Mass  and  Differentiation  matrices  for 
%2  -  Dimension . 

%Written  by  Benjamin  Davis  Created:  July  2014  Modified  Jan  2015 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  % 

function  [M,D]  =  mass_diff2D  (L,dL,w) 
n  =  length (L ( : ,1) ) ; 

M  =  zeros (n,n) ; 

D  =  zeros (n,n) ; 
for  i  =  l:n 


for  j  =  1 :  n 

M(  i  ,  j )  =  ( ( L  ( i  , : )  ,*L(j  ,:))*w(l  ,:)  ')  ; 
D(  i  ,  j )  =  L(i  ,:)  *dL(j  ,:) 

end 

end 

end 
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%Written  by  Benjamin  Davis 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  23  Jan  2015 

%Synopsis  :  Building  Matrix  Terms  in  order  to  solve  the  equations  for 
%the  2D  Acoustic  Wave  Equation . 

%Output :  Matrix  elements  used  to  solve  the  RHS  in  the  RK54  scheme. 

% - % 

function  [MAT, MAT_ksi , MAT_eta ,  Je  ,A,B,C,D,E,F,G,PI]  =  MatrixTerms2D  (Ne, 
D_ksi ,  D_eta  ,  y_eta  ,  y_ksi  ,  x_eta  ,  x_ksi  ,  jac  ,M) 
for  k  =  l:Ne 

Je { k }  =  diag (jac (k , : ) ) ; 

Jeinv  { k }  =  diag  ( 1  ./ jac  (k  ,:))  ; 
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A  =  D_ksi  '  X-  (  diag  ( y_eta  (k  , 
B  =  D_eta  '  *  (  diag  ( y_ksi  (k  , 
C  =  D_ksi  '  *  (  diag  ( y_eta  (k  , 
D=  D_eta  '  *  (  diag  ( y_ksi  (k  , 
E  =  D_ksi  '  *  (  diag  (  x_eta  (k  , 
F  =  D_eta  '  *  (  diag  (  x_ksi  (k  , 
G  =  D_ksi  '  x  (  diag  (  x_eta  (k  , 
H=  D_eta  '  x  (  diag  (  x_ksi  (k  , 


)  )  x  Jeinv  {k }  xMx  diag  (  y_eta  (k  , 
)  )  x  Jeinv  {k  }  xMx  diag  (  y_eta  (k  , 
)  )  x  Jeinv  {k  }xMx  diag  (  y_ksi  (k  , 
)  )  x  Jeinv  {k  }  xMx  diag  (  y_ksi  (k  , 
)  )  x  Jeinv  {k  }  xMx  diag  (  x_eta  (k  , 
)  )  x  Jeinv  {k  }  xMx  diag  (  x_eta  (k  , 
)  )  x  Jeinv  { k  }  xMx  diag  (  x_ksi  (k  , 
)  )  x  Jeinv  {k }  xMx  diag  (  x_ksi  (k  , 


MAT{k[  =  A-B-C+DfE-F-G+H; 

MAT_ksi  { k  j  =  A-B+E-F; 

MAT_eta  { k  j  =  -C+D-G+H;%%%G+D-G+H 
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end 


end 


)  ) )  xD_ksi ; 
)  ) )  xD_ksi ; 
)  )  )  xD_eta  ; 
)  )  )  xD_eta  ; 
)  ) )  xD_ksi ; 
)  ) )  xD_ksi ; 
)  )  )  xD_eta  ; 
)  )  )  xD_eta  ; 
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o/  o/ 

/o - /o 

%Function  to  build  Dx  and  Dy  used  to  solve  the  RHS  in  the  RK54  scheme. 
%Written  by  Benjamin  Davis 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  23  Jan  2015 

% - % 

function  [Dx,Dy]  =  DxDy(Ne,  ksi_x  ,  ksi_y  ,  eta_x  ,  eta_y  ,  D_ksi ,  D_eta ) 
for  k  =  l:Ne 


Dxjkj  =  diag  (  ksi_x  (k  , :  )  )  xD_ksi  +  diag  (  eta_x  (k  , : )  )  xD_eta  ; 
Dyjkj  =  diag  (  ksi_y  (k  , :  )  )  xD_ksi  +  diag  (  eta_y  (k  , : )  )  xD_eta  ; 

end 

end 
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0/  o/ 

/o - /o 

%Function  for  building  the  Surface  Jacobians  for  the  faces  on  a  2D  grid 
%  Written  by  Benjamin  Davis 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  23  Jan  2015 

% - % 

function  [  Sjl  ,  Sj 2  ,  Sj 3  ,  S j 4  ]  =  SurfaceJac2D  (Ne,Ll , L2 , L3  ,L4 ,  x_ksi  ,  x_eta  , 
y_ksi , y_eta ) 
for  k  =  l:Ne 


%Building  the  surface  jacobians  for  the  faces 
S j 1 _3 { k }  =  sqrt ( x_ksi (k , : )  . A2  +  y_ksi (k , : )  . A2) ; 

Sj 2_4  { k  }  =  sqrt  (  x_eta  (k  , : )  .  A2  +  y_eta  (k  , : )  .  A2)  ; 
%Isolate  Surface  Jacobians  individually  for  four  faces 
Sjl  { k (  =  diag(Ll*Sjl_3  { k }  ')  ; 

S j 3  { k }  =  diag(L3*Sjl_3  { k }  ')  ; 

S j 2  { k }  =  diag(L2*Sj2_4  { k }  ')  ; 

S j 4  { k }  =  diag(L4*Sj2_4  jk)  ')  ; 

end 

end 
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o/  o/ 

/o - /o 

%Function  to  build  the  element  normals  for  the  fours  faces  on  a  2D  grid 
%Written  by  Benjamin  Davis 

%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

%  23  Jan  2015 

% - % 


function  [  nx_l  ,  nx_2  ,  nx_3  ,  nx_4  ,  ny_l  ,  ny_2  ,  ny_3  ,  ny_4  ]  =  ElementNormals2D  (Ne 
,L1  ,L2,L3,L4,  x_ksi  ,  x_eta  ,y_ksi  ,y_eta  ,  Sjl  ,  S j 2  ,  Sj 3  ,  S j 4  ) 
for  k  =  l:Ne 

%Compute  the  Normals  per  element 
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nx_l 

ik) 

=  diag  (LI  *y_ksi  (k  , : ) 

•  /  diag  (  Sj  1  { k  j )  )  ; 

ny_l 

!  k } 

=  -  diag  ( LI  *  x_ksi  (k  , : ) 

'  ■/  diag  (  S j  1  { k } ) ) ; 

nx_2 

!  k  j 

=  diag  (L2*y_eta  (k  , : ) 

./ diag  (  Sj  2  { k } ) )  ; 

ny_2 

{ k } 

=  -diag (L2* x_eta (k , : ) 

' ./ diag  (  Sj 2  { k } )  )  ; 

nx_3 

!  k } 

=  -diag  (L3*y_ksi  (k  , : ) 

' ./ diag  (  Sj 3  { k } )  )  ; 

ny_3 

!  k  j 

=  diag  ( L3  *  x_ksi  (k  , : ) 

./ diag  (  Sj  3  { k  } )  )  ; 

nx_4 

!  k } 

=  -diag  (L4*y_eta (k , : ) 

'  •/  diag  (  S j 4  { k  } )  )  ; 

ny_4 { k } 

=  diag  ( L4* x_eta (k , : ) 

■/  diag  (  S j  4  { k } ) ) ; 

end 

end 
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%Function  provided  to  Professor  F.X.  Giraldo  '  s  MA4245  class. 
%Modified  by  Asst.  Professor  Jeremy  Kozdon  and  Benjamin  Davis 
%This  function  rotates  the  grid  and  plots  it  . 

%Written  by  F.X.  Giraldo  on  4/2008 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% 

%  INPUT  LIST :  coord  are  the  coordinates 
%  intma  is  the  connectivity  list 

%  npoin  are  the  number  of  global  points 

%  nelem  are  the  number  of  elements 

%  ngl  is  the  number  of  interpolation  points  in  an  element 

%  (polynomial  order  +  1)\ 

%  plot_grid  is  a  switch  to  either  plot  or  not  plot 

%  grid_rotation_angle  is  the  grid  rotation  in  degrees 

%  OUTPUT  LIST : 

%  coord  are  the  new  rotated  coordinates:  x=coord(:,l) 

%  and  y=coord(:,2) 

O/  O/ 

/o - /o 
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function  [  coord_rotated  ]  =  rotate_grid_v2  ( coord  ,  intma  ,  npoin  ,  nelem  ,  ngl , 
plot_grid  ,  grid_rotation_angle ) 
nop=ngl - 1 ; 

coord_rotated  =  zeros  (npoin  ,2  )  ; 

%Rotate  Grid 

alpha  =  grid_rotation_angle*pi/180; 
for  i=l:npoin 

xn=cos ( alpha ) *coord ( i , 1 )  -  sin ( alpha ) *  coord ( i ,2 ) ; 

yn=sin ( alpha ) *coord ( i , 1 )  +  cos ( alpha ) *  coord ( i ,2 ) ; 
coord_rotated(i ,1) =xn ; 
coord_rotated(i ,2) =yn ; 
end 

%Plot  Grid 
if  (plot_grid  ==  1) 
x=zeros (5,1) ; 
y=zeros (5,1)  ; 
figure  ; 
hold  on; 
for  e  =  l:nelem 

for  j  =l:ngl  -1 

for  i=l:ngl-l 

il  =intma  ( e  ,  i  ,  j  )  ; 
i2=intma  ( e  ,  i  +1 ,  j  )  ; 
i3=intma (e,  i  +1, j  +1) ; 
i4=intma  (e  ,  i  ,  j  +1) ; 

x(  1 )  =coord_rotated  ( il  ,1)  ;  y ( 1 )  =coord_rotated  ( il  ,2)  ; 
x(2)=coord_rotated ( i2 ,1)  ;  y(2)=coord_rotated ( i2  ,2)  ; 
x(3)=coord_rotated ( i3 ,1)  ;  y(3)=coord_rotated ( i3  ,2)  ; 
x(4)=coord_rotated(i4 ,1)  ;  y(4)=coord_rotated(i4  ,2)  ; 
x(5)=coord_rotated  ( il  ,1)  ;  y(5)=coord_rotated  ( il  ,2)  ; 
plot_handle=plot (x,y, ' -r ' ) ; 
set  ( plot_handle  ,  '  Line  Width  '  ,1.5)  ; 

end 

end 

i  1  =intma  ( e ,  1  , 1 )  ; 
i2=intma  (e  ,  ngl  ,  1 )  ; 
i3=intma  (e  ,  ngl  ,  ngl )  ; 
i4=intma(e,l  ,ngl)  ; 

x( 1 )=coord ( il ,1) ;  y ( 1 ) =coord ( il ,2) ; 
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x(2) =coord ( i2  ,1)  ;  y (2 ) =coord ( i2  ,2)  ; 
x(3) =coord ( i3  ,1)  ;  y (3) =coord ( i3 ,2)  ; 
x(4)=coord(i4 ,1)  ;  y (4) =coord ( i4 ,2 ) ; 
x(5) =coord ( il  ,1) ;  y (5 ) =coord ( il  ,2)  ; 
plot_handle  =  plot (x,y,  '  -b  '  ) ; 
set  ( plot_handle  ,  '  LineWidth  ,2)  ; 

end 

title_text  =[  ' Rotated  Grid  Plot  For:  Ne  =  '  num2str(nelem)  N  = 
num2str(nop)  ]; 

title  ([  title_text  ]  ,  'FontSize  '  ,18) ; 

xlabel(  'X'  ,  'FontSize  '  ,18)  ; 

ylabel(  'Y'  ,  'FontSize  '  ,18)  ; 

axis  image 

end 
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%Function  for  building  the  Lagrange  Polynomials . 

%Written  by  Benjamin  Davis  in  MA4245  %Created  :  July  2014 
%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - %  % 

function  [L,dL]  =  lagrange_basis  ( x ,  z  ) 

%Nth  order  interpolation 
n  =  length ( x ) ; 

%Length  of  the  equally  spaced  grid  for  k  =  1:50 
h  =  length ( z )  ; 

%Initialize  the  Lagrange  Matrix 
L  =  ones (n,h) ; 
dL  =  zeros (n,h) ; 

%Computation  for  Lagrange  Matrix 
for  k  =  l:h 

for  i  =  l:n 

for  j  =  l:n 

dl  =  1; 
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if  j  ~=  i  %  If  j  does  not  equal  i 

%Equation  for  the  Lagrange  Polynomial 
L(i,k)  =  (z(k)-x(j)) ./(x(i)-x(j))  *  L(i,k); 
for  1  =  l:n 

if  (1  ~=  i)  &&  (1  -=  j ) 

dl  =  dl*(z(k)-x(l)) ./(x(i)-x(l)); 

end 

end 

dL(i,k)  =  dL ( i , k )  +  dl/(x(i)-x(j )); 

end 

end 

end 

end 

end 
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% - % 

%Code  given  to  F.X.  Giraldo  MA4245  Class  July  2014 
%Used  by  Ben  Davis  . 

%This  code  computes  the  Legendre -Gauss - Lobatto  points  and  weights 
%which  are  the  roots  of  the  Lobatto  Polynomials . 

%Written  by  F.X.  Giraldo  on  4/2000 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 


function  [xgl,wgl]  =  legendre_gauss_lobatto  (P) 
p=P - 1 ; 

ph=floor  (  (p  +  l)/2  ); 
for  i=l:ph 

x=cos(  (2*  i -1 )  *pi /( 2*p  +  l)  ); 
for  k  =  l:20 

[L0,L0_1  ,L0_2]  =  legendre_poly  (p,x) ; 

%Get  new  Newton  Iteration 

dx=-(l-xA2)*L0_l/(  -2*x*L0_l  +  ( 1 -xA2) *L0_2 ) ; 
x=x+dx ; 

if  (abs(dx)  <  1.0e-20) 
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break 

end 

end 

X§1  (p+2-i  )=x; 

wgl(p+2-i )  =  2/(p*(p  +  l)  *L0 A2) ; 

end 

“/(Check  for  Zero  Root 
if  (p+1  ~=  2*ph) 
x  =  0; 

[L0,L0_1  ,L0_2]  =  legendre_poly  (p,x) ; 
xgl  (ph  +  l)=x; 

wgl(ph  +  l)  =  2/(p*(p  +  l)*L0A2) ; 

end 

%Find  remainder  of  roots  via  symmetry 
for  i=l:ph 

xgl(i)=-xgl(p+2-i); 
wgl(  i  )=+wgl(p+2-i )  ; 

end 
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% - % 

%Code  given  to  F.X.  Giraldo  MA4245  Class  July  2014 
%Used  by  Ben  Davis 

%This  code  computes  the  Legendre  Polynomials  and  its  1st  and  2nd 
%derivatives 

% 

%Written  by  F.X.  Giraldo  on  4/2000 


%  Department  of  Applied  Mathematics 

%  Naval  Postgraduate  School 

%  Monterey,  CA  93943-5216 

% - % 


function  [L0 ,  L0_1 ,  L0_2 ]  =  legendre_poly  (p,x) 
L1=0;L1_1=0;L1_2=0; 

L0  =  1;L0_1  =0;L0_2=0; 
for  i=l:p 

L2=L1 ; L2_1=L1_1 ; L2_2=L1_2 ; 

L1=L0 ; L1_1=L0_1 ; L1_2=L0_2 ; 
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a  =  (2*i  -1)  /  i ; 

b=(i  -1)  /  i ; 

LO=a*x*Ll  -  b*L2; 

LO_l=a*(Ll  +  x*Ll_l)  -  b*L2_l; 
L0_2=a*(2*Ll_l  +  x*Ll_2)  -  b*L2_2; 

end 
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